Hi all,

I am solving a homogenization problem on a composite consisting of 2 different materials. I have defined the composite with 2 subdomains on one mesh. I just wanted to know how I can imply continuity of stress and displacement at the interface? does it count as a boundary condition in fenics?

If the subdomains are topologically-connected in the mesh, then solving for displacement in a CG function space defined on this mesh will strongly enforce displacement continuity at the interface. If you just vary the material properties on the two subdomains within a standard Galerkin weak form, that will also weakly enforce traction compatibility.

To see this, suppose the subdomains are \Omega_1 and \Omega_2, with \Omega = \Omega_1\cup\Omega_2. Galerkinâs method poses the weak problem

where \pmb{\sigma}_i is the Cauchy stress in \Omega_i. Integrating by parts on each \Omega_i, one gets

where \Gamma = \overline{\Omega}_1\cap\overline{\Omega}_2, i.e., the shared boundary or interface (and Iâve assumed that strongly-enforced Dirichlet and/or periodic BCs are enforced on any external boundaries), and \mathbf{n}_i is the outward-facing normal to \Omega_i. Because \mathbf{n}_1 = -\mathbf{n}_2 on \Gamma, the EulerâLagrange equation corresponding to the \Gamma integral is

which is the stress continuity that you want physically, for equilibrium, while the volume integrals correspond to the strong form of -\nabla\cdot\pmb{\sigma}_i = \mathbf{f} on each subdomain.

Thanks a lot for your reply. It was very helpful.

So, as you said, I just used CG function space but I get the following result which is shown in the picture below. I can tell that stress contour is discontinuous. Am I right?

It depends on what scalar quantity youâre plotting as âstressâ. Equilibrium doesnât imply that all components of the stress tensor are continuous; you only expect to have (\pmb{\sigma}_1 - \pmb{\sigma}_2)\cdot\mathbf{n} = \mathbf{0}.

Thatâs right. But this is the contour for y-axis normal stress.

I just simplified my code to check out if I made a mistake in my variational formulation, but I got the same discontinuity in my result. Here is the simplified code that I just wrote and the final result that I got.

Thanks in advance.

for name in dir():

if not name.startswith(â_â):

del globals()[name]

from dolfin import *

import numpy as np

import matplotlib.pyplot as pl

class Top(SubDomain):

def inside(self, x, on_boundary):

return near(x[1], 1.0)

top = Top()

class Left1(SubDomain):

def inside(self, x, on_boundary):

return near(x[0], 0.0) and between(x[1], (0.5, 1))

leftup=Left1()

class Right1(SubDomain):

def inside(self, x, on_boundary):

return near(x[0], 1.0) and between(x[1], (0.5, 1))

rightup = Right1()

class Interface(SubDomain):

def inside(self, x, on_boundary):

return near(x[1], 0.5)

interface = Interface()

class Left2(SubDomain):

def inside(self, x, on_boundary):

return near(x[0], 0.0) and between(x[1], (0.0, 0.5))

leftdown = Left2()

class Right2(SubDomain):

def inside(self, x, on_boundary):

return near(x[0], 1.0) and between(x[1], (0.0, 0.5))

rightdown = Right2()

class Bottom(SubDomain):

def inside(self, x, on_boundary):

return near(x[1], 0.0)

bottom = Bottom()

#Create meshes

n = 50

mesh = UnitSquareMesh(n, n)

class Solidup(SubDomain):

def inside(self, x, on_boundary):

return x[1] >= 0.5

solidup= Solidup()

class Soliddown(SubDomain):

def inside(self, x, on_boundary):

return x[1] <= 0.5

soliddown = Soliddown()

regions = MeshFunction(âsize_tâ, mesh, mesh.topology().dim())

regions.set_all(0)

solidup.mark(regions, 1)

soliddown.mark(regions, 2)

boundaries = MeshFunction(âsize_tâ, mesh, 1)

boundaries.set_all(0)

top.mark(boundaries, 1)

leftup.mark(boundaries, 2)

rightup.mark(boundaries, 3)

interface.mark(boundaries, 4)

leftdown.mark(boundaries, 5)

rightdown.mark(boundaries, 6)

bottom.mark(boundaries, 7)

#define a submesh, composed of region 0 and 1

solidup = MeshFunction(âsize_tâ, mesh, mesh.topology().dim())

solidup.set_all(0)

solidup.array()[regions.array() == 1] = 1

submeshsolidup = SubMesh(mesh, solidup, 1)

soliddown = MeshFunction(âsize_tâ, mesh, mesh.topology().dim())

soliddown.set_all(0)

soliddown.array()[regions.array() == 2] = 2

submeshsoliddown = SubMesh(mesh, soliddown, 2)

#define new meshfunctions on this submesh with the same values as their original mesh

solidup_regions = MeshFunction(âsize_tâ, submeshsolidup, submeshsolidup.topology().dim())

solidup_regions.set_all(0)

solidup_boundaries = MeshFunction(âsize_tâ, submeshsolidup, 1)

solidup_boundaries.set_all(0)

soliddown_regions = MeshFunction(âsize_tâ, submeshsoliddown, submeshsoliddown.topology().dim())

soliddown_regions.set_all(0)

soliddown_boundaries = MeshFunction(âsize_tâ, submeshsoliddown, 1)

soliddown_boundaries.set_all(0)

W = VectorFunctionSpace(mesh,âCGâ, 2) # Electroelasticity

w = Function(W)

uz = TestFunction(W)

dx = Measure(âdxâ, domain=mesh, subdomain_data=regions)

bcs = [DirichletBC(W, Constant((0.,0.)), bottom)

,DirichletBC(W, Constant((0.,0.4)), top)]

E = Constant(1e5)

nu = Constant(0.3)

model = âplane_stressâ

Landa1=100

Landa2=150

mu = E/2/(1+nu)

lmbda = Enu/(1+nu)/(1-2nu)

if model == âplane_stressâ:

lmbda = 2mulmbda/(lmbda+2*mu)

def sigma(v):

return lmbdatr(grad(v))grad(v)Identity(2) + 2.0mu

Wint =inner(Landa1*sigma(uz), grad(w))

dx(1)+ inner(Landa2sigma(uz), grad(w))*dx(2)

solve(Wint==0, w, bcs)

#Plot sigma and u

plt.figure()

p3 = plot(w[1],title=âksiâ)

plt.colorbar(p3)

plt.show()

plt.figure()

ps = plot(grad(w)[1,1],title=ây-axis Normal stress [MPa]â)

plt.colorbar(ps)

plt.show()

First, I believe thereâs a typo in your formulation, since, in the line `Wint = ...`

, it would be more standard for `uz`

to be the `Function`

and `w`

to be the `TestFunction`

, and for the nonlinear solve to have `uz`

as its unknown solution. Second, if `w`

is interpreted as a displacement, `grad(w)[1,1]`

is not the normal stress in the y-direction, it is the extensional strain. (In general, strain would be the symmetrized gradient of displacement, but symmetrization doesnât change diagonal components.) You do not expect strain to be continuous. The stress tensor would be `Landa1*sigma(uz)`

or `Landa2*sigma(uz)`

, depending on which subdomain youâre in (assuming the code is changed so that `uz`

is the displacement, rather than the test function).

Just to add to what David has already said, if your function is in `CG2`

(piecewise quadratic approximation), then itâs gradient should be in `DG1`

(piecewise linear discontinuous). This is what you should see in your results. I donât know right off the bat, but `plot`

may/may not handle `DG1`

functions. The best way to visualize is by writing them using write_checkpoint to a `xdmf`

file and open in Paraview.

Also, others would be able to help you better, if you enclose your code withing triple fences `````

such that it is reproducible.

Youâre totally right. I had some typos in that code, sorry for that. Now, I totally get it. just a further question if you donât mind. As I am new to fenics, I have no idea that how I will be able to plot sigma(uz) in each subdomain individually. I mean how to separate sigma(uz) in each subdomain to multiply to different Landa.

Thanks for your time.

Thanks for your reply. I got the xdmf file and the result which is shown in paraview is same as the one that I got in by using plot. As david said I have to plot uz in in each subdomain to use the right coefficient for that. I just wanted to ask you if you know how to plot sigma(uz) in each subdomain individually. I mean how to separate sigma(uz) in each subdomain to multiply to different Landa.

So, the way the existing `plot`

function works is that most UFL objects you pass as the first argument are projected (in the sense of L^2) onto a CG1 (linear, continuous) finite element space for visualization. (If you pass a CG1 function, it is plotted directly, or, if you pass a DG0 function, that can also be plotted directly.) You can mimic this behavior for a piecewise function that is `f1`

on subdomain 1 and `f2`

on subdomain 2 by writing a custom projection function:

```
def split_project(f1,f2,V):
u = TrialFunction(V)
v = TestFunction(V)
uh = Function(V)
solve(u*v*dx == f1*v*dx(1) + f2*v*dx(2),uh)
return uh
```

The return value of this will be a `Function`

in the space `V`

that is the L^2 projection of the piecewise function defined by `f1`

and `f2`

. If `V`

is a CG1 or DG0 space, that can be plotted with `plot`

directly. (Otherwise, it will be projected onto CG1 to plot.) As @bhaveshshrimali pointed out, you can get a higher-fidelity visualization by choosing `V`

to be a higher-degree DG space, but you would need to write it out to a file and plot with Paraview.

An extra note: L^2-projecting a discontinuous function onto a CG space will give oscillatory results; you see this in your original plot of `grad(w)[1,1]`

(which was projected onto CG1, as explained above), where there is some over- and under-shoot near the discontinuity. One way around such oscillatory behavior is to use a lumped-mass L^2 projection, as demonstrated in my answer here; that is also more computationally-efficient than solving a linear system for the consistent L^2 projection.

How about using a `DG0`

function `landa`

for that takes values `Landa1`

and `Landa2`

as you wish; then using `landa*sigma(uz)`

as your function to be (local-) projected onto a `DG1`

space and finally written to `Paraview`

using `write_checkpoint`

., which in your case should roughly translate to

```
from dolfin import *
metadata={"quadrature_degree": 2}
def local_project(v,V):
dv = TrialFunction(V)
v_ = TestFunction(V)
a_proj = inner(dv,v_)*dx(metadata=metadata)
b_proj = inner(v,v_)*dx(metadata=metadata)
solver = LocalSolver(a_proj,b_proj)
solver.factorize()
u = Function(V)
solver.solve_local_rhs(u)
return u
V1 = TensorFunctionSpace(mesh, "DG", 1) #for sigma
Vdg = FunctionSpace(mesh, "DG", 0) # for landa
landa = Function(Vdg)
x = Vdg.tabulate_dof_coordinates()
for i in range(x.shape[0]):
if regions.array()[i] == 1:
landa.vector().vec().setValueLocal(i, 100) # `Landa1`
else:
landa.vector().vec().setValueLocal(i, 150) # `Landa2`
sigma_vis = local_project(landa*sigma(uz), V1)
with XDMFFile(MPI.comm_world, "sigma.xdmf") as fil:
fil.write_checkpoint(sigma_vis, "sigma", 0)
```

Thanks a lot.That works perfectly.

Great.That works perfectly.

Please, you could post the complete code as bhaveshshrimali commented, all of it encapsulated with ``` !

And have you remade considering the application of a load on the upper domain?

Dear kamensky,

Could you point me to a book or bibliographic reference where you got this approach from the two-domain garlekin method?

The basic idea is probably as old as the structural analysis that modern finite element methods historically grew out of, since structural analysts are happy to use adjacent elements of different stiffness. I based the particular notation I used in my post on this paper about fluidâstructure interaction.