Surface function evaluated at unexpected points

Hi,
I am solving a Poisson equation with a non-linear surface density on an internal surface.
The code runs but while doing few crosschecks I noticed that the non-linear function is evaluated also for points outside the internal surface.
Does anyone understand this behaviour?
Here is the implementation of the problem

V = dol.FunctionSpace(
    mesh_dic['mesh'], 
    "P", 1,  #linear Lagrange elements
    )
# trial and test functions
u = dol.Function(V)
v = dol.TestFunction(V)
du = dol.TrialFunction(V)
# measures 
dx = mesh_dic['volume_measure']
ds = mesh_dic['boundary_measure']
dS = mesh_dic['internal_boundary_measure']
# quadratic form and its jacobian
a = dol.inner(dol.grad(u), dol.grad(v)) * dx
da = dol.inner(dol.grad(du), dol.grad(v)) * dx

nl = dol.avg(v) * Rho(u) * dS(3)
dnl = dol.avg(v) * dol.avg(du) * Rho_Der(u) * dS(3)
bc = []
bc.append(
    dol.DirichletBC(V,
                    dol.Constant(0),
                    mesh_dic['boundary_marker'],
                    1))
bc.append(
        dol.DirichletBC(V,
                    dol.Constant(1),
                    mesh_dic['boundary_marker'],
                    2)
)
    # solve non-linear problem
functional = a - nl
jacobian = da - dnl
problem = dol.NonlinearVariationalProblem(functional, u, bc, jacobian)
solver = dol.NonlinearVariationalSolver(problem)
solver.solve()

and here the classes I used to implement the non linear functions

class Rho(dol.UserExpression):
    def __init__(self, phi):
        self.phi = phi
        dol.UserExpression.__init__(self)
    def eval(self, value, x):
        print(x[2]) # I added this to check where the function is evaluated
        mu = self.phi(x) 
        value[0] = - mu**2 * np.sign(mu) 
    def value_shape(self):
        return ()
class Rho_Der(dol.UserExpression):
    def __init__(self, phi):
        self.phi = phi
        dol.UserExpression.__init__(self)
    def eval(self, value, x):
        mu = self.phi(x) 
        value[0] = - 2. * abs(mu) 
    def value_shape(self):
        return ()

Thanks a lot.
Iacopo

The question is not just academic because for more complicated geometries the function gets evaluated outside the mesh resulting in an error.

If the positions are only slightly off the surface, it may just be some sort of floating-point precision issue, which you can resolve by including

self.phi.set_allow_extrapolation(True)

in the constructor.

Also, you might consider implementing eval_cell(self,value,x,cell) instead of eval(self, value, x), if you only plan to use Rho in form assembly (where you would call self.phi.eval_cell instead of using the overloaded __call__ method of self.phi to eval it). This is generally more efficient in that it doesn’t require the computational geometry step of figuring out what cell the point x is in.

1 Like

The points are off by quite significant amounts (around the size of the finite elements) it is not machine precision.

I already tried to allow for extrapolation and, indeedd, the code works but I am not sure that it is doing the right thing.

The function is actually defined on facets, is there an equivalent of eval_cell for facets?

The following code snippet prints out \sqrt{2}, as expected from interior facet integration, so I think you can just use eval_cell for facets as well as cells:

from dolfin import *
mesh = UnitSquareMesh(1,1)
class Test(UserExpression):
    def eval_cell(self,value,x,cell):
        value[0] = 1.0
print(assemble(Test()*dS(domain=mesh)))

Thanks a lot!
But this still doesn’t make me feel totally sure about what the code is doing.

Hi, does anybody (perhaps @dokken) have an idea of what is going on?
I still do not understand why the function is evaluated outside the surface.
Best. Iacopo

You need to post a minimal example that reproduces the issue (that me or others can run on our system and get the same issues as you have encountered).
The code above does not include a mesh, as it is dependent on the mesh.dic, and should be written as a single code block, that can be copy-pasted.

import os
import numpy as np
import matplotlib.pyplot as plt
import dolfin as dol
import meshio

def convert_msh_to_XDMF(filename, dim, dim_keep = None):
    '''reads a .msh file and writes mesh, subdomains, and boundary subdomains 
    to XDMF files'''
    if filename.split(sep = '.')[-1] != 'msh':
        raise TypeError('.msh file required')
    fn = '.'.join(filename.split(sep ='.')[:-1])
    if dim_keep is None:
        dim_keep = dim
    msh = meshio.read(filename)
    if dim == 1:
        volume_cell_type = 'line'
        surface_cell_type = 'vertex'
    elif dim == 2:
        volume_cell_type = 'triangle'
        surface_cell_type = 'line'
    elif dim == 3:
        volume_cell_type = 'tetra'
        surface_cell_type = 'triangle'
    else:
        raise ValueError('dim can only be 1,2, or 3')
    
    volume_cells = []
    surface_cells = []
    for cell in msh.cells:
        if cell.type == volume_cell_type:
            if len(volume_cells) == 0:
                volume_cells = cell.data
            else:
                volume_cells = np.vstack([volume_cells, cell.data])

        if cell.type == surface_cell_type:
            if len(surface_cells) == 0:
                surface_cells = cell.data
            else:
                surface_cells = np.vstack([surface_cells, cell.data])
        
    volume_data = None
    surface_data = None
    for key in msh.cell_data_dict['gmsh:physical'].keys():
        if key == volume_cell_type:
            volume_data = msh.cell_data_dict['gmsh:physical'][key]
        elif key == surface_cell_type:
            surface_data = msh.cell_data_dict['gmsh:physical'][key]
    volume_mesh = meshio.Mesh(points = msh.points[:,:dim_keep],
                                cells = [(volume_cell_type, volume_cells)],
                                cell_data = {'marker' : [volume_data]})
    meshio.write(fn + '_{}D_mesh.xdmf'.format(dim), volume_mesh)
    if not (surface_data is None):
        surface_mesh = meshio.Mesh(points = msh.points[:,:dim_keep], 
                                    cells = [(surface_cell_type, surface_cells)],
                                    cell_data = {'marker' : [surface_data]})
        meshio.write(fn + '_{}D_mesh.xdmf'.format(dim-1), surface_mesh)
    

def load_XDMF_files(filename, dim):
    '''create mesh from XDMF files'''
    mesh = dol.Mesh()
    fn = '.'.join(filename.split(sep ='.')[:-1])
    with dol.XDMFFile(fn + '_{}D_mesh.xdmf'.format(dim)) as infile:
        infile.read(mesh)
    mvc = dol.MeshValueCollection('size_t', mesh = mesh, dim = dim) 
    with dol.XDMFFile(fn + '_{}D_mesh.xdmf'.format(dim)) as infile:
        infile.read(mvc, 'marker')
    volume_marker = dol.cpp.mesh.MeshFunctionSizet(mesh, mvc)
    del(mvc)
    dx = dol.Measure('dx', domain = mesh, subdomain_data = volume_marker)

    if os.path.exists(fn + '_{}D_mesh.xdmf'.format(dim-1)):
        mvc = dol.MeshValueCollection('size_t', mesh = mesh, dim = dim -1) 
        with dol.XDMFFile(fn +  '_{}D_mesh.xdmf'.format(dim-1)) as infile:
            infile.read(mvc, 'marker')
        boundary_marker = dol.cpp.mesh.MeshFunctionSizet(mesh, mvc)
        del(mvc)
        ds = dol.Measure('ds', domain = mesh, subdomain_data = boundary_marker)
        dS = dol.Measure('dS', domain = mesh, subdomain_data = boundary_marker)
    else:
        warn('{}D mesh not found'.format(dim-1))
        boundary_marker = None
        ds = dol.Measure('ds', domain = mesh)
        dS = dol.Measure('dS', domain = mesh)
    return {'mesh' : mesh, 
            'volume_measure' : dx,
            'boundary_measure' : ds,
            'internal_boundary_measure' : dS,
            'volume_marker' : volume_marker,
            'boundary_marker' : boundary_marker}

class Rho(dol.UserExpression):
    def __init__(self, phi):
        self.phi = phi
        dol.UserExpression.__init__(self)
    def eval(self, value, x):
        print(x[2]) # I added this to check where the function is evaluated
        mu = self.phi(x) 
        value[0] = - mu**2 * np.sign(mu) 
    def value_shape(self):
        return ()
class Rho_Der(dol.UserExpression):
    def __init__(self, phi):
        self.phi = phi
        dol.UserExpression.__init__(self)
    def eval(self, value, x):
        mu = self.phi(x) 
        value[0] = - 2. * abs(mu) 
    def value_shape(self):
        return ()

#Actual program starts here

convert_msh_to_XDMF( 'structure.msh', dim = 3)
mesh_dic = load_XDMF_files('structure.msh', dim = 3)

V = dol.FunctionSpace(
    mesh_dic['mesh'], 
    "P", 1,  #linear Lagrange elements
    )
# trial and test functions
u = dol.Function(V)
v = dol.TestFunction(V)
du = dol.TrialFunction(V)
# measures 
dx = mesh_dic['volume_measure']
ds = mesh_dic['boundary_measure']
dS = mesh_dic['internal_boundary_measure']
# quadratic form and its jacobian
a = dol.inner(dol.grad(u), dol.grad(v)) * dx
da = dol.inner(dol.grad(du), dol.grad(v)) * dx

nl = dol.avg(v) * Rho(u) * dS(3)
dnl = dol.avg(v) * dol.avg(du) * Rho_Der(u) * dS(3)
bc = []
bc.append(
    dol.DirichletBC(V,
                    dol.Constant(0),
                    mesh_dic['boundary_marker'],
                    1))
bc.append(
        dol.DirichletBC(V,
                    dol.Constant(1),
                    mesh_dic['boundary_marker'],
                    2)
)
    # solve non-linear problem
functional = a - nl
jacobian = da - dnl
problem = dol.NonlinearVariationalProblem(functional, u, bc, jacobian)
solver = dol.NonlinearVariationalSolver(problem)
solver.solve()

Here is the code but to run it you will need the .msh file.
How can I attach it?
Best.
Iacopo

Please make sure your problem is a minimal working example, that only has enough code to illustrate your issue, not the whole application. I would either suggest you attach the geo file, or try to solve a similar problem with a built in mesh, and check if you observe similar behavior.

Hi, here is the .geo file.
I am sorry but the issue is happening on internal boundaries and I don’t know how to generate a mesh with internal boundaries without using an external program.
The biggest part of what I posted is just to import the externally generated mesh in Fenics.
It comes from a script that I wrote following your instructions and always worked fine.
Iacopo

// translation vectors
a1_x = 48.0;
a1_y = 0.0;
a2_x = 24.0;
a2_y = 41.569219381653056;

// hole 1
d1_x = 36.0;
d1_y = 20.784609690826528;
r1 = 5.0;

// layer 1
t_1 = 50;

// layer 2
t_2 = 1;

// layer 3
t_3 = 3.5;

// layer 4
t_4 = 20;

SetFactory("OpenCASCADE");

// create bottom face
Point(1) = {0,0,0};
Point(2) = {a1_x,a1_y,0};
Point(4) = {a2_x, a2_y, 0};
Point(3) = {a1_x + a2_x, a1_y + a2_y, 0};
Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {3,4};
Line(4) = {4,1};
Curve Loop(5) = {1,2,3,4};
Plane Surface(1) = {5};

// first box
Extrude {0,0,t_1} {Surface{1};}

//cylinder
Cylinder(2) = {d1_x, d1_y,t_1,0,0,t_2, r1};

// second box
Translate {0, 0, t_1 +t_2} { Duplicata{Surface{1};}}
Extrude {0,0,t_3} {Surface{10};}

// third box
Translate {0, 0, t_1 +t_2+t_3} { Duplicata{Surface{1};}}
Extrude {0,0,t_4} {Surface{16};}

// intersections
BooleanFragments{Volume{1,2};Delete;}{}
BooleanFragments{Volume{2,3};Delete;}{}
BooleanFragments{Volume{3,4};Delete;}{}

// physical volumes
Physical Volume('substrate', 1) = {1};
Physical Volume('holes', 2) = {2};
Physical Volume('bottom_layer', 3) = {3};
Physical Volume('top_layer', 4) = {4};

// physical surfaces
Physical Surface('bottom_gate', 1) = {26};
Physical Surface('patterned_gate', 2) = {27,7,32};
Physical Surface('graphene', 3) = {33};
Physical Surface('top_gate', 4) = {38};

Mesh.MinimumCirclePoints = 20;
Field[1] = Box;
Field[1].VIn = 2;
Field[1].VOut = 50;
Field[1].XMin = 0.0;
Field[1].XMax = a1_x +a2_x;
Field[1].YMin = 0.;
Field[1].YMax = a1_y +a2_y;
Field[1].ZMin = 50;
Field[1].ZMax = 54.5;
Background Field = 1;

My point with the previous posts, is that you can illustrate the behavior you are seeing with alot less code, i.e. consider:

from dolfin import *

mesh = UnitCubeMesh(1,1,10)
V = FunctionSpace(mesh, "CG", 1)

mf = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)

class Internal(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[2], 0.5)

internal_boundary = Internal()
internal_boundary.mark(mf, 3)
File("test_mf.pvd") << mf


class coord(UserExpression):
    def __init__(self):
        UserExpression.__init__(self)
    def eval(self, value, x):
        print(x) # I added this to check where the function is evaluated
        value[0] = x[2] 
    def value_shape(self):
        return ()


dS_c = Measure("dS", domain=mesh, subdomain_data=mf, metadata={"quadrature_degree": 4})
print(assemble(coord()*dS_c(3)))
print(assemble(1*dS_c(3)))

This returns a set of points (that are not on the surface), but the correct integral value:

0.49999999999999956
0.999999999999999

This can be explained by the fact that even if you are integrating over a facet, the UserExpression is interpolated into an appropriate function space behind the scenes. The function space lives on the cells, and as a consequence, every cell that is connected to the interface will have an evaluation at them, to get the correct expression.
Note that this is less expensive than manually interpolating the function into a space, as it would then do it for all cells (as you can verify with interpolate(coord(), V)).

Hi, thanks for your reply.
I understand that the function gets interpolated by using evaluations outside the integration surface (even thought it is not obvious to me that this can be an optimal strategy) but this should happen under the hood and I would not have noticed it if it wasn’t giving an error.
Using the code (and mesh) I posted I get the function evaluated at points that are not only outside the surface but also outside the function domain, resulting in an error.
The issue can be solved by allowing extrapolation, I guess it is fine but I was scared that the result may be affected in some way.

Note, moreover that if I try to reproduce the behaviour by using only an integration (I should have tried this before posting) I don’t encounter any issue (see attached code). The problem appears only when the nonlinear solver is called.

import os
import numpy as np
import matplotlib.pyplot as plt
import dolfin as dol
import meshio

def convert_msh_to_XDMF(filename, dim, dim_keep = None):
    '''reads a .msh file and writes mesh, subdomains, and boundary subdomains 
    to XDMF files'''
    if filename.split(sep = '.')[-1] != 'msh':
        raise TypeError('.msh file required')
    fn = '.'.join(filename.split(sep ='.')[:-1])
    if dim_keep is None:
        dim_keep = dim
    msh = meshio.read(filename)
    if dim == 1:
        volume_cell_type = 'line'
        surface_cell_type = 'vertex'
    elif dim == 2:
        volume_cell_type = 'triangle'
        surface_cell_type = 'line'
    elif dim == 3:
        volume_cell_type = 'tetra'
        surface_cell_type = 'triangle'
    else:
        raise ValueError('dim can only be 1,2, or 3')
    
    volume_cells = []
    surface_cells = []
    for cell in msh.cells:
        if cell.type == volume_cell_type:
            if len(volume_cells) == 0:
                volume_cells = cell.data
            else:
                volume_cells = np.vstack([volume_cells, cell.data])

        if cell.type == surface_cell_type:
            if len(surface_cells) == 0:
                surface_cells = cell.data
            else:
                surface_cells = np.vstack([surface_cells, cell.data])
        
    volume_data = None
    surface_data = None
    for key in msh.cell_data_dict['gmsh:physical'].keys():
        if key == volume_cell_type:
            volume_data = msh.cell_data_dict['gmsh:physical'][key]
        elif key == surface_cell_type:
            surface_data = msh.cell_data_dict['gmsh:physical'][key]
    volume_mesh = meshio.Mesh(points = msh.points[:,:dim_keep],
                                cells = [(volume_cell_type, volume_cells)],
                                cell_data = {'marker' : [volume_data]})
    meshio.write(fn + '_{}D_mesh.xdmf'.format(dim), volume_mesh)
    if not (surface_data is None):
        surface_mesh = meshio.Mesh(points = msh.points[:,:dim_keep], 
                                    cells = [(surface_cell_type, surface_cells)],
                                    cell_data = {'marker' : [surface_data]})
        meshio.write(fn + '_{}D_mesh.xdmf'.format(dim-1), surface_mesh)
    

def load_XDMF_files(filename, dim):
    '''create mesh from XDMF files'''
    mesh = dol.Mesh()
    fn = '.'.join(filename.split(sep ='.')[:-1])
    with dol.XDMFFile(fn + '_{}D_mesh.xdmf'.format(dim)) as infile:
        infile.read(mesh)
    mvc = dol.MeshValueCollection('size_t', mesh = mesh, dim = dim) 
    with dol.XDMFFile(fn + '_{}D_mesh.xdmf'.format(dim)) as infile:
        infile.read(mvc, 'marker')
    volume_marker = dol.cpp.mesh.MeshFunctionSizet(mesh, mvc)
    del(mvc)
    dx = dol.Measure('dx', domain = mesh, subdomain_data = volume_marker)

    if os.path.exists(fn + '_{}D_mesh.xdmf'.format(dim-1)):
        mvc = dol.MeshValueCollection('size_t', mesh = mesh, dim = dim -1) 
        with dol.XDMFFile(fn +  '_{}D_mesh.xdmf'.format(dim-1)) as infile:
            infile.read(mvc, 'marker')
        boundary_marker = dol.cpp.mesh.MeshFunctionSizet(mesh, mvc)
        del(mvc)
        ds = dol.Measure('ds', domain = mesh, subdomain_data = boundary_marker)
        dS = dol.Measure('dS', domain = mesh, subdomain_data = boundary_marker)
    else:
        warn('{}D mesh not found'.format(dim-1))
        boundary_marker = None
        ds = dol.Measure('ds', domain = mesh)
        dS = dol.Measure('dS', domain = mesh)
    return {'mesh' : mesh, 
            'volume_measure' : dx,
            'boundary_measure' : ds,
            'internal_boundary_measure' : dS,
            'volume_marker' : volume_marker,
            'boundary_marker' : boundary_marker}

class coord(dol.UserExpression):
    def __init__(self):
        dol.UserExpression.__init__(self)
    def eval(self, value, x):
        if x[2] != 54.5:
            print(x) # I added this to check where the function is evaluated
        value[0] = x[2] 
    def value_shape(self):
        return ()

#Actual program starts here

convert_msh_to_XDMF( 'structure.msh', dim = 3)
mesh_dic = load_XDMF_files('structure.msh', dim = 3)
dS = mesh_dic['internal_boundary_measure']

print(dol.assemble(coord()*dS(3)))
print(dol.assemble(1*dS(3))*54.5)