Problem with MeshFunction

Question

Hi folks,

Just a quick question about MeshFunction. I am trying to have three different sub-meshes (mesha, meshs and meshc in the code)and function spaces (V_a, V_s, V_c).

I use the following code to generate these. But when I want to print out ca(Lxa), I get the following error. I absolutely have no idea why this is happening! Lxa is inside the domain, not outside.

Any help is deeply appreciated.

Thanks,
Mo

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

nxa = 20
nxs = 5
nxc = 20

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 a(SubDomain):
def inside(self, x, on_boundary):
return True if 0 <= x[0] <= (Lxa) else False

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

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

subdomain0 = a()
subdomain0.mark(subdomains, 0)
subdomain1 = s()
subdomain1.mark(subdomains, 1)
subdomain2 = c()
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)

c = Function(V)
ca = Function(V_a)
cc = Function(V_c)

print ca(0)
print ca(Lxa)
print cc(Lxa+Lxs)
print cc(Lx)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

And the 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

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

Answer 1

Hi,

Maybe this is caused by the point Lxa being not inside the mesha. You can try to modify the subdomain definition of s and c if the codes work after you replace ca(Lxa) by cc(Lxa).

Comments

Hi Mirok,

Thank you for your reply. But why this is not happening at point Lxa + Lxs, when I want to print cs(Lxa + Lxs) and cc(Lxa + Lxs)? Isn't it strange to you?

Is there any way we can modify this such that Lxa belongs to both mesh0 and mesh1? So far, near, between and introducing a tolerance did not solve the issue!

Thanks,
Mo

---

Hi, near and between were part of your inside methods, but as I said, these are not used
in function evaluation. There, mesh.intersected_cell is called but I am not sure if and how you
could introduce some tolerance to calculations of that function. Anyways, is allowing extrapolation a problem in your application?

---

Hi Mirok,

My problem is not evaluating a function. The issue here is that I do not think I am using a correct mesh and I think this is the reason of getting different solutions when I change the number of elements. Lxa should belong to both mesha and meshc, otherwise, I will not get a correct solution. Don't you agree?

May I report a bug? Do we consider this a bug? The reason I am saying this is that the issue does not appear at point Lxa + Lxs. Another strange thing is that it works for timesn=90! It really surprises me!

Thanks,
Mo

---

Hi, I get error running your code even with timesn=90. I think the differences you are seeing between Lxa and Lxa+Lxs are due to finite precision but I see no harm in reporting this behavior.

---

Thank you Mirok. I do appreciate your time and your continued help. I will report this.

Best,
Mo