Mixed Dimensional Branch: lagrange multiplier

variational formulation
1

Hi,
I am trying to solve a linear homogenization problem as described in here but using the total displacement as the unknown.

TLDR; I am solving a simple linear elasticity problem, looking to constrain the displacement on the right (u_R) and left (u_L) edges of a simple unit square as:

\hspace{15em} u_R - u_L = \widetilde{u} (prescribed)

I am thinking of using a lagrange multiplier (\lambda \in \mathbb{R}^2 (\Gamma)) defined on the left and right edge to impose the above constrain, and tweaking the bilinear and linear forms respectively as:

\hspace{1em} a((u, \lambda), (\delta u, \delta\lambda)) = {\large\int}_\Omega \sigma(u):\varepsilon(\delta u)dx + {\large \int}_\Gamma\ \ \delta\lambda\cdot (\delta u_R - \delta u_L) d\Gamma + {\large \int}_\Gamma\ \ \lambda \cdot (\delta u_R - \delta u_L) d\Gamma

\hspace{5em} L(\delta u, \delta\lambda) = {\large\int}_\Gamma \delta\lambda\cdot\widetilde{u} d\Gamma

I am using @cdaversin 's mixed-dimensional branch.

Here is a simple MWE that fails:

from dolfin import *
mesh = UnitSquareMesh(100, 100)


left = CompiledSubDomain("near(x[0], 0.) && on_boundary")
top = CompiledSubDomain("near(x[1], 1.) && on_boundary")
bottom = CompiledSubDomain("near(x[1], 0.) && on_boundary")
right = CompiledSubDomain("near(x[0], 1.) && on_boundary")

origin = CompiledSubDomain("near(x[0], 0.) && near(x[1], 0.) && on_boundary")

facetMarks= MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
funcMarks = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)

E, nu = Constant(2.), Constant(0.25)
mu = E/2./(1+nu)
lmbda = 2*mu*nu/(1-2*nu)
eps = lambda x: sym(grad(x))
sigma = lambda x: lmbda*tr(eps(x)) * Identity(len(x)) + 2*mu*eps(x)

facetMarks.set_all(0)
funcMarks.set_all(0)

bottom.mark(facetMarks, 1)
top.mark(facetMarks, 2)
left.mark(facetMarks, 3)
right.mark(facetMarks, 4)

bottom.mark(funcMarks, 1)
top.mark(funcMarks, 1)
left.mark(funcMarks, 2)
right.mark(funcMarks, 2)

mesh_B = MeshView.create(facetMarks, 1)
mesh_T = MeshView.create(facetMarks, 2)
mesh_L = MeshView.create(facetMarks, 3)
mesh_R = MeshView.create(facetMarks, 4)

mesh_LR = MeshView.create(funcMarks, 2)
mesh_TB = MeshView.create(funcMarks, 1)

Wu = VectorElement("CG", mesh.ufl_cell(), 2)
Wlam_LR = VectorElement("R", mesh_LR.ufl_cell(), 0)
Wlam_TB = VectorElement("R", mesh_TB.ufl_cell(), 0)

V_u = FunctionSpace(mesh, Wu)
V_lam_LR = FunctionSpace(mesh_LR, Wlam_LR)
V_lam_TB = FunctionSpace(mesh_TB, Wlam_TB)

V = MixedFunctionSpace(V_u, V_lam_LR, V_lam_TB)

u, lam_R, lam_T = TrialFunctions(V)
uh, lam_Rh, lam_Th = TestFunctions(V)

# volume measures
ds_R = Measure("dx", domain=mesh_R)
ds_L = Measure("dx", domain=mesh_L)
ds_T = Measure("dx", domain=mesh_T)
ds_B = Measure("dx", domain=mesh_B)

# dirichelt BC's
bcs = [DirichletBC(V.sub_space(0), Constant((0,0)), origin, method="pointwise"),
       DirichletBC(V.sub_space(0).sub(1), Constant(0.0), CompiledSubDomain("near(x[0], 1.) && near(x[1], 0.0) && on_boundary"), method="pointwise"),
       DirichletBC(V.sub_space(0).sub(0), Constant(0.0), CompiledSubDomain("near(x[0], 0.) && near(x[1], 1.0) && on_boundary"), method="pointwise") ]

a = inner(sigma(u), eps(uh))*dx
a += inner(lam_Rh, u)* ds_R - inner(lam_Rh, u)*ds_L + inner(lam_R, uh)* ds_R - inner(lam_R, uh)* ds_L
a += inner(lam_Th, u)* ds_T - inner(lam_Th, u)* ds_B + inner(lam_T, uh)* ds_T - inner(lam_T, uh)* ds_B

L = inner(lam_Rh, Constant(("0.1", "0"))) * ds_R + inner(lam_Th, Constant(("0.", "-0.1")))* ds_T

usol = Function(V, name="u")

solve(a == L, usol, bcs, solver_parameters={"linear_solver": "lu"})

I get the following error:

PetscDolfinErrorHandler: line '503', function 'MatSetValues_SeqAIJ', file '/tmp/petsc-src/src/mat/impls/aij/seq/aij.c',
                       : error code '63' (Argument out of range), message follows:
------------------------------------------------------------------------------
New nonzero at (80782,0) caused a malloc
Use MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE) to turn off this check
2

@cdaversin Have you had a chance to look at this?

3

Hi,
Thanks for your interest in the mixed dimensional branch.
I am not a 24/7 support so no, I haven’t taken the time to look at this over the weekend.
I will try to have a look ASAP, and I’ll keep you posted.

2 Likes
4

Hi,

The problem here is that the Measure you are using in the form integrals are not based on the same domain (=mesh) as the trial and/or test functions involved :
For instance, in inner(lam_Rh, u)* ds_R : lam_Rh is in V_lam_LR based on mesh_LR, u is in V_u based on mesh but the measure ds_R is defined from domain = mesh_R.
Your measure should be based on mesh_LR instead, possibly with some subdomain_data.

But I guess you’ve already figured that out, considering your other post here.

1 Like
5

Thanks

I tried tweaking the variational formulation to have Measures that make more sense in the current setting. Whenever you get time, you may look if you are able to reproduce the error in the linked post.

1 Like