setting condition for mesh refinement

Question

Hi,

I am trying to set a condition for where to refine the mesh and I'm getting the following error message:

if phi0(p) < 0.05 and phi0(p) > -0.05:
File "C:\FEniCS\lib\site-packages\dolfin\functions\function.py", line 342, in
__call__
raise TypeError, "expected scalar arguments for the coordinates"
TypeError: expected scalar arguments for the coordinates

I have obtained the function phi0 by solving a PDE, and I would like to refine the mesh in all cells where this function attains values between -0.05 and 0.05. The corresponding piece of code is:

cell_markers = CellFunction("bool", mesh)
cell_markers.set_all(False)
for cell in cells(mesh):
    p = cell.midpoint()
    if phi0(p) < 0.05 and phi0(p) > -0.05:
        cell_markers[cell] = True
mesh = refine(mesh, cell_markers)

Should I be addressing the cells in a different way? I could use a subdomain definition to cover the range for phi0 but the problem is it would overlap with my defined subdomains. Is there a way out?

Many thanks,

Answer 1

Your code is good. This is bug in DOLFIN which is already resolved in development version. The workaround is to call phi0([p.x(), p.y(), p.z()]) instead of phi0(p).

But note that if you would need to optimize it for efficiency you should use phi0.eval(value, x, cell) because you already have an information on which cell you want to evaluate phi0 and you waste this information and let geometry algorithm to search again for cell. This note does not apply iff phi0 is defined on different mesh.

Comments

Sure. Thanks very much.

On another note (for mesh refinement), is there a way I could set how small I would like the refined mesh elements to be?

I understand the method in place is edge bisection, and the way I see the refined mesh, it has divided the original elements into two which is far too coarse for my computations. If I repeat the

mesh = refine(mesh, cell_markers)

command, would that do the job of re-dividing?

Many thanks,

---

I'm using the latest version of fenics from the latest Docker image. I am trying to do something similar to the original question. I tried using the eval method per this answer; but the arguments seem to be wrong. I get this error:

NotImplementedError: Wrong number or type of arguments for overloaded function 'Function_eval'.
Possible C/C++ prototypes are:
dolfin::Function::eval(dolfin::Array< double > &,dolfin::Array< double > const &) const
dolfin::Function::eval(dolfin::Array< double > &,dolfin::Array< double > const &,dolfin::Cell const &,ufc::cell const &) const

So my short question is: What was assumed to be the type of phi0 in the original question? What has changed in the past four years that change the answer?

Longer question: Should I be doing this differently if my solution Function is from a mixed space (made with FunctionSpace using a mixed element)?

I somewhat new to fenics and I have trouble navigating the documentation. I spent a bit of time trying to answer my question by reviewing the code, but I'm not getting anywhere. Thanks for the help.

Edit/Update: If I use eval_cell instead of eval, I instead run into the error here: https://answers.launchpad.net/fenics/+question/207055

Edit/Update 2: Now I'm making some progress. I was failing to use Jan's advice about eval. First off, in the suggested call, phi0.eval(value, x, cell), "x" must not be the same type as "p" from the question (a fenics.Point), but rather it should be x=[p.x(), p.y(), p.z()], which was mentioned as a "workaround" for a bug that was supposedly fixed in a 2013 development version. Furthermore, regarding the error message from my original comment, the "eval" method does not appear to have the correct arguments, but "eval_cell" seems to work. So now I have a running line of code which is "w.eval_cell(values, numpy.array([p.x(), p.y(), p.z()]), cell)".