Hello everyone,
I am solving an interface problem with Nitsche’s method where I want to impose stress continuity, and displacement continuity at the interface of two elastic bodies. I’ve found an excellent example here (Stress continuity using mixed-dimensional branch) that works very well with Lagrange multipliers. Only I’ve chosen Nitsche’s method because I don’t want to add extra variables to penalise the problem.
I’m not familiar with Nitsche’s method and how it works. I can’t always understand how the interface is managed. Should I create the interface before making the mesh or should I split the mesh a bit like in the attached code. Is the best solution to solve for a global displacement u and therefore use the jump and avg functions in FEniCS?
As mentioned above, I started from an example already shared here which I tried to adapt using the Nitsche method. The code is below. I’m getting an error that I can’t fix.
I need your expertise to understand this method and how to implement it with FEniCS. Sorry, I’m new to FEniCSx and as I have a deadline, I wanted to test it with FEniCS first.
Thanks in advance.
from dolfin import *
class Dirichlet(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 2) \
or near(x[1], 0) \
or near(x[1], 1)
class Neumann(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0)
class Interface(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1)
# Mesh generation
n = 20
mesh = RectangleMesh(Point(0,0,0),Point(2,1,0),2*n,n)
# Domain markers
domain = MeshFunction("size_t", mesh, mesh.geometry().dim(), 0)
for c in cells(mesh):
domain[c] = c.midpoint().x() > 1.0
# Interface markers
facets = MeshFunction("size_t", mesh, mesh.geometry().dim()-1, 0)
Interface().mark(facets,1)
# Meshes
mesh_L = MeshView.create(domain, 0)
mesh_R = MeshView.create(domain, 1)
mesh_I = MeshView.create(facets, 1)
# Function spaces
FE = VectorElement("CG", mesh_I.ufl_cell(), 1, dim=2) # Lagrange multiplier
VE = VectorElement("CG", mesh.ufl_cell(), 1, dim=2) # Displacement
V_L = FunctionSpace(mesh_L,VE) # Displacement on left domain
V_R = FunctionSpace(mesh_R,VE) # Displacement on right domain
V_I_disp = FunctionSpace(mesh_I,FE) # Lagrange multiplier - displacement continuity
V_I_stress = FunctionSpace(mesh_I,FE) # Lagrange multiplier -stress continuity
V_I = FunctionSpace(mesh_I,FE)
W = MixedFunctionSpace(V_L, V_R)
# Function and test functions
(ul, ur) = TrialFunctions(W)
(wl, wr) = TestFunctions(W)
# Markers for Dirichlet and Neuman bcs
ff_L = MeshFunction("size_t", mesh_L, mesh_L.geometry().dim()-1, 0)
Dirichlet().mark(ff_L, 1)
Neumann().mark(ff_L, 2)
Interface().mark(ff_L, 3)
ff_R = MeshFunction("size_t", mesh_R, mesh_R.geometry().dim()-1, 0)
Dirichlet().mark(ff_R, 1)
Interface().mark(ff_R, 3)
# Boundary conditions
bc_L = DirichletBC(W.sub_space(0), Constant((0.0,0.0)), ff_L, 1)
bc_R = DirichletBC(W.sub_space(1), Constant((0,0)), ff_R, 1)
bcs = [bc_L, bc_R]
# Elasticity parameters
rho = Constant(2200)
lmbda = Constant(7.428e+04)
mu = Constant(7.428e+04)
# Neumann condition
g = Constant((1e+5, 0))
# Measures
dx_L = Measure("dx", domain=W.sub_space(0).mesh())
ds_L = Measure("ds", domain=W.sub_space(0).mesh(), subdomain_data=ff_L)
dx_R = Measure("dx", domain=W.sub_space(1).mesh())
ds_R = Measure("ds", domain=W.sub_space(1).mesh(), subdomain_data=ff_R)
#ds_I = Measure("dx", domain=W.sub_space(2).mesh())
# Définition de l'interface
class Interface(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1)
# Marquage de l'interface sur le maillage
facet_marker = MeshFunction("size_t", mesh_L, mesh_L.topology().dim() - 1)
facet_marker.set_all(0)
interface = Interface()
interface.mark(facet_marker, 4)
#ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
# Stress tensor
def sigma(r):
return 2.0*mu*sym(grad(r)) + lmbda*tr(sym(grad(r)))*Identity(len(r))
# Variational problem
a_left = inner(sigma(ul),sym(grad(wl)))*dx_L + inner(sigma(ul)*FacetNormal(mesh_L), wl)*ds_L(2)
L_left = inner(g, wl)*ds_L(2)
a_right = inner(sigma(ur),sym(grad(wr)))*dx_R #- inner(sigma(ul)*FacetNormal(mesh_R), wr)*ds_R(3)
L_right = inner(Constant((0,0)), wr)*dx_R
nL = Constant((1,0))
nR = Constant((-1,0))
# definition fo the Nitsche
Thet = -1;
gamma = 10
lam = 0.5
h= CellDiameter(mesh_L)
avg_u = lam*sigma(ul)+ (1-lam)*sigma(ur)
avg_w = lam*sigma(wl)+ (1-lam)*sigma(wr)
jum_u = ul-ur
jum_w = wl-wr
a = a_left + a_right - inner(avg_u*nL, jum_w)*ds_L(3) + Thet*inner(avg_w*nL, jum_u)*ds_L(3) + 100* inner(jum_u,jum_w)*ds_L(3)
L = L_left + L_right
# Output files
pvda_file = File("output1/ul.pvd")
pvdb_file = File("output1/ur.pvd")
#pvdc_file = File("output/uc.pvd")
# Solve
u_mixed = Function(W)
solve(a == L, u_mixed, bcs, solver_parameters={"linear_solver":"direct"})
ul, ur= u_mixed.sub(0), u_mixed.sub(1)
# Save solution
pvda_file << ul
pvdb_file << ur
#pvdc_file << uc
I have the following error,
*** -------------------------------------------------------------------------
*** Error: Unable to successfully call PETSc function 'MatSetValuesLocal'.
*** Reason: PETSc error code is: 63 (Argument out of range).
*** Where: This error was encountered inside /build/dolfin-pnfrym/dolfin-2019.2.0~git20201207.b495043/dolfin/la/PETScMatrix.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: unknown
*** -------------------------------------------------------------------------
Thanks in Advance
