Using MeshFunction instead of FacetFunction to define ds measure

I would like to replace FacetFunction by MeshFunction in order to comply with the new version of FEniCS. The old version of my code is

w = 50
L = 100.
mesh = RectangleMesh(Point(0., 0.), Point(L, w), 5, 5, "crossed")

def eps(v):
    return sym(grad(v))

mu = Constant(1.0)
lmbda = Constant(0.5)

def sigma(v):
    return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)

class Right(SubDomain):
    def inside(self, x, on_boundary):
	    return near(x[0], L)

exterior_facet_domains = FacetFunction("size_t", mesh)
exterior_facet_domains.set_all(0)
Right().mark(exterior_facet_domains, 3)
ds = Measure('ds', domain=mesh, subdomain_data=exterior_facet_domains)
s = Constant((1, 0.0))

V = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
l = inner(s, u_)*ds(3)


def left(x, on_boundary):
    return near(x[0], 0.)

def bottom(x, on_boundary):
    return near(x[1], 0.)

bcx = DirichletBC(V.sub(0), Constant((0.)), left)
bcy = DirichletBC(V.sub(1), Constant((0.)), bottom)
bc = [bcx, bcy]

u = Function(V)
solve(a == l, u, bc)

Now when I replace

exterior_facet_domains = FacetFunction("size_t", mesh) 

by

exterior_facet_domains = MeshFunction("double", mesh, mesh.topology().dim()-1)

I get the following error:

TypeError: set_exterior_facet_domains(): incompatible function arguments. The following argument types are supported:
1. (self: dolfin.cpp.fem.Form, arg0: dolfin.cpp.mesh.MeshFunctionSizet) -> None

Any help would be very welcome!

Change the "double" argument to "size_t" and it should work

that was easy, thanks!

Hello Dokken,

I have a very similar problem but can’t figure out why. I also get “set_exterior_facet_domains(): incompatible function arguments” as my failure. This comes in combination with “Memory access error (memory dumped)”

from dolfin import *

if not has_petsc():
    print("DOLFIN must be compiled at least with PETSc 3.6 to run this demo.")
    exit(0)
    
parameters["linear_algebra_backend"] = "PETSc"

mesh = UnitSquareMesh(10,10)

V = VectorFunctionSpace(mesh, "Lagrange", 1)
u, du, v 	= Function(V), TrialFunction(V), TestFunction(V)
x = SpatialCoordinate(mesh)

class left(SubDomain):
    def inside(self, x, on_boundary):
    	return x[0] < 0.001

def eps(u):
    return sym(grad(u))
		
constraint_l 	= Expression(("z_min-x[0]", "r_min-x[1]"), r_min=0, z_min=0, degree=2)
constraint_u 	= Expression(("z_max-x[0]", "r_max-x[1]"), r_max=0.9, z_max=1.4, degree=2)
u_min 		= interpolate(constraint_l, V)
u_max 		= interpolate(constraint_u, V)

l_mf = MeshFunction("size_t", mesh, mesh.topology().dim() - 1) 
l_mf.set_all(0)
left().mark(l_mf, 1) 
ds = Measure('ds', domain = mesh, subdomain_data = left)

class LEFT(UserExpression):
	def _init_(self, mesh, **kwargs):
		self.mesh=mesh
		super()._init_(**kwargs)	
	def eval(self, values, x):
		values[0] = 1
		values[1] = 1
	def value_shape(self):
		return (2,)

f= LEFT(mesh)

E, nu 	= 8* 10**6, 0.4
mu 	= Constant(E/(2.0*(1.0+nu)))
lmbda 	= Constant(E*nu/((1.0+nu)*(1.0-2.0*nu)))

I 	= Identity(2)
F 	= I + eps(u)
C 	= F.T*F
Ic	= tr(C)
J 	= det(F)
psi	= (mu/2)*(Ic-3)-mu*ln(J)+(lmbda/2)*(ln(J))**2

elastic_energy 	= psi*dx - dot(f, u)*ds
grad_elastic_energy	= derivative(elastic_energy, u, v)
H_elastic_energy 	= derivative(grad_elastic_energy, u, du)

    def __init__(self):
        OptimisationProblem.__init__(self)


    def f(self, x):
        u.vector()[:] = x
        return assemble(elastic_energy)

    def F(self, b, x):
        u.vector()[:] = x
        assemble(grad_elastic_energy, tensor=b)

    def J(self, A, x):					
        u.vector()[:] = x
        assemble(H_elastic_energy, tensor=A)
       
solver = PETScTAOSolver()

solver.parameters["method"] = "bqpip"	
solver.parameters["monitor_convergence"]	= True
solver.parameters["gradient_absolute_tol"]	= 1.0e-1 #1.0e-08 <--default
solver.parameters["gradient_relative_tol"]	= 1.0e-06 #1.0e-08 <--default
solver.parameters["report"] 			= False
solver.parameters["maximum_iterations"]	= 10000

parameters.parse()

solver.solve(BucklingProblem(), u.vector(), u_min.vector(), u_max.vector())

I think I got that, ‘ds’ is predefined as a measurement over the whole surface of my mesh. And as long as I do not redefine it and function f is integrated over the whole surface everything works fine.

But redefining the measurement ‘ds’ on the subdomain ‘left’ causes the error.

Please help :slight_smile:

The subdomain_data should be l_mf, not left.

You should probably use ds(1) and not ds when defining the variational form if you only want to integrate over the left facets

1 Like