Problem with MeshFunction

Hi folks,

I am trying to have three 1D sub-meshes (mesha, meshs and meshc in the code) and function spaces (V_a, V_s, V_c) using MeshFunction.

I use the following code to generate these. The issue is that point x=Lxa is supposed to blong to both subdomain0 and subdomain1. Similarly, x=Lxa+Lxs should belong to both subdomain1 and subdomain2. Although x=Lxa+Lxs belongs to both subdomain1 and subdomain2, but point x=Lxa only blongs to subdomain1.

When I want to print out ca(Lxa), I get the following error.

Traceback (most recent call last):
File "test_mesh.py", line 71, in
print ca(Lxa)
File "/usr/lib/python2.7/dist-packages/dolfin/functions/function.py", line 382, in call
self.eval(values, x)
RuntimeError:

*** -------------------------------------------------------------------------

*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at

*** https://answers.launchpad.net/dolfin

*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.

*** -------------------------------------------------------------------------
*** Error: Unable to evaluate function at point.
*** Reason: The point is not inside the domain. Consider setting "allow_extrapolation" to allow extrapolation.
*** Where: This error was encountered inside Function.cpp.
*** Process: 0

The reason I am reporting this as a bug is that there is no problem at point x=Lxa+Lxs, while at x=Lxa, the issue appears.

I have also posted a qustion on FEniCS Q & A which you may follow the whole story by refering to the following link:

http://fenicsproject.org/qa/1479/problem-with-meshfunction?show=1486#a1486

My code is:

import time
from dolfin import *
import numpy as np
import scipy.io as sciio
from math import tanh, sinh, cosh, log
from numpy import linalg as LA
import scipy.io
set_log_active(False)

from sys import exit

Lxa = 100.0E-6                
Lxs = 25.0E-6               
Lxc = 100.0E-6      

Lx = Lxa + Lxs + Lxc 

timesn = 10

nxa = 4*timesn             
nxs = 1*timesn             
nxc = 4*timesn            

nx = nxa + nxs + nxc         

dxa = Lxa/nxa
dxs = Lxs/nxs
dxc = Lxc/nxc

print 'dxa ===', dxa
print 'dxs ===', dxs
print 'dxc ===', dxc

mesh = IntervalMesh(nx,0,Lx)

subdomains = MeshFunction('uint', mesh, 0)

class ano(SubDomain):
    def inside(self, x, on_boundary):
        return True if between(x[0], (0, Lxa)) else False

class sep(SubDomain):
    def inside(self, x, on_boundary):
        return True if between(x[0], (Lxa, Lxa+Lxs)) else False

class cat(SubDomain):
    def inside(self, x, on_boundary):
        return True if between(x[0], (Lxa+Lxs, Lx)) else False

# Mark subdomains with numbers
subdomain0 = ano()
subdomain0.mark(subdomains, 0)
subdomain1 = sep()
subdomain1.mark(subdomains, 1)
subdomain2 = cat()
subdomain2.mark(subdomains, 2)

mesha = SubMesh(mesh, subdomain0)    
meshs = SubMesh(mesh, subdomain1)    
meshc = SubMesh(mesh, subdomain2)   


V = FunctionSpace(mesh, 'CG', 1)         
V_a = FunctionSpace(mesha, 'CG', 1)         
V_s = FunctionSpace(meshs, 'CG', 1)         
V_c = FunctionSpace(meshc, 'CG', 1)     

mesha_coor = mesha.coordinates()   
meshc_coor = meshc.coordinates()     

ca = Function(V_a)
cs = Function(V_s)
cc = Function(V_c)

print  'Mesh coordinates in ano:', mesha_coor
print  'Mesh coordinates in cat:', meshc_coor

print 'in ano at foil:', ca(0)
print 'in ano at separator:', ca(Lxa)  
print 'in sep at separator:', cs(Lxa)    
print 'in sep at separator:', cs(Lxa+Lxs)
print 'in cat at separator:', cc(Lxa+Lxs)
print 'in cat at foil:', cc(Lx)

Any help is deeply appreciated.

Thanks,

Mo

Comments (4)