Define integral on part of boundary in energy formulation

Question

Hi,

I am solving a variationnal problem in non linear elasticity (Neo-Hookean model) involving an external force F acting on a part of the boundary (indexed by 2 inMeshFunction exterior_facets_meshfunction) :

mesh = Mesh("mymesh.xml")
exterior_facets_meshfunction = MeshFunction("size_t", mesh, 'mymeshfunction.xml")
V = VectorFunctionSpace(mesh, "Lagrange", 1)
u  = Function(V)
ds = Measure("ds")[exterior_facets_meshfunction]
load = Expression(("t"),t=F)

I then define the invariants, material parameters and energy density psi.

I tried two ways to define the total energy (taking into account the effect of the external load on boundary) :

W = psi*dx - load*u[1]*ds(2)(mesh)

and :

W = psi*dx - load*u[1]*ds(2)

Does someone know the difference between the two syntaxes ?

I have the impression that the first formulation does the integration over the whole mesh, so I would rather use the second one, but it raises the following error message when I try to define the first variation of energy using derivative :

ufl.log.UFLException: An Integral without a Domain is now illegal.

Is there another way to restrict integration on a part of the boundary ?

I am using version 1.7.0dev

Many thanks in advance !

Claire

Answer 1

Hi, here's an explanation of what happens

from dolfin import *

mesh = UnitSquareMesh(10, 10)
facet_f = FacetFunction('size_t', mesh, 0)
CompiledSubDomain('near(x[0], 0)').mark(facet_f, 2)

ds = Measure("ds")[facet_f]

V = FunctionSpace(mesh, 'CG', 1)
c = Expression('1', domain=mesh)

# Are the measures the same?
print assemble(c*ds(2))
print assemble(c*ds(2)(mesh))

# No! So what happens 
print ds(2).subdomain_id()          # 2 as expected
print ds(2)(mesh).subdomain_id()    # everywhere, why?

# Because ds(2)(mesh) is a call to ds(2).__call__ and that results in
print ds(2).__call__.__doc__ 

Comments

Hi, many thanks for your answer.

So I understand the two measures are different, and that the second one does the integration over the whole mesh, so that's not the one I want to use.

Nevertheless I obtain an error when I use the first one in the following code (classical non linear hyper elasticity) :

from dolfin import *

mesh = UnitSquareMesh(10, 10)
facet_f = FacetFunction('size_t', mesh, 0)
CompiledSubDomain('near(x[0], 0)').mark(facet_f, 2)

ds = Measure("ds")[facet_f]

V = VectorFunctionSpace(mesh, 'CG', 1)

mu, Kpar = 1., 10.
bcs = []

# define residual 
def F(u, lmbda, v):
    I = Identity(2)
    FF = I + grad(u)
    C = FF.T*FF
    Ic = tr(C)
    JJ  = det(FF)
    psi = (mu/2.*(Ic-2 - 2*ln(JJ))+Kpar/2.*(ln(JJ))**2)*dx - lmbda*u[1]*ds(2)
    return derivative(psi, u, v)

# define variational problem
u0 = Function(V)

dE = F(u0, 0., TestFunction(V))
ddE=derivative(dE,u0,TrialFunction(V))
problem=NonlinearVariationalProblem(dE,u0,bcs,J=ddE)
solver=NonlinearVariationalSolver(problem)

The error appears when dE is defined :

ufl.log.UFLException: An Integral without a Domain is now illegal.

How can I solve this issue ? Many thanks in advance !!

---

Where is load defined?

---

Sorry I mixed up load and lmbda, replace by :

psi = (mu/2.*(Ic-2 - 2*ln(JJ))+Kpar/2.*(ln(JJ))**2)*dx - lmbda*u[1]*ds(2)

(I made the correction in my previous comment...)

---

Try with ds = Measure("ds", domain=mesh, subdomain_data=facet_f)

---

It works fine ! Many thanks again for your help !