Local damaged elastic material

Local damaged elastic material#

Initialize#

For parallel testing#

Create a server controlling n>1 MPI engines (consider OpenMPI as MPI backend)

#nopy
import ipyparallel as ipp
nw=2
cluster = ipp.Cluster(engines="mpi", n=nw)
rc = cluster.start_and_connect_sync()
view=rc[:]

/dolfinx-env/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm
[stderr:0] 
/dolfinx-env/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm
[stderr:1] 
/dolfinx-env/lib/python3.12/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

DOLFINx version: 0.11.0.post0 based on GIT commit: fefdb2201b80a8f59527de2d461b9056906661d8

Import TwoScale library#

%%px --local
from twoscale import core
from twoscale import linear
from twoscale import util

Use Twoscale dolfinx implementation
Use Twoscale dolfinx implementation

[stdout:0] 
Use Twoscale dolfinx implementation
Use Twoscale dolfinx implementation
[stdout:1] 
Use Twoscale dolfinx implementation
Use Twoscale dolfinx implementation

Basic test case#

2026-09-23 15:48:34.563 (   1.009s) [    7F8185B21140]vtkXOpenGLRenderWindow.:1450  WARN| bad X server connection. DISPLAY=
../../_images/aaf1301aa8559729a22a77a89e19beae2a75628b6350b967e96733765c98fc86.png 
[stderr:0] 2026-09-23 15:48:35.011 (   1.606s) [    7F1801CA4140]vtkXOpenGLRenderWindow.:1450  WARN| bad X server connection. DISPLAY=

[output:0]
../../_images/5cfad0df04129be130992672a9644c5b081d9b7e713dcd559bbca517e6c0b7fb.png

Choose a mesh refinement localization criterion#

finest mesh in an oblique band in the middle

%%px --local
delta=3*L/(ref*4*s)
def crit1(coords,level):
    return np.logical_and(coords[0]<L*(coords[1]/B+s)/(2*s)+ref*delta/3,coords[0]>L*(coords[1]/B+s)/(2*s)-ref*delta/3)

Choose macro node to enrich#

select central region around x=L/2 +/- L/8

%%px --local
def enriched(coords):
    return np.logical_and(coords[0]>=3*L/8, coords[0]<=5*L/8)

Scale jump construction#

%%px
j=core.topDown(domain,enriched,crit1,ref)

%%px --local
fdomain=j.getFineMesh
cdomain=j.getCoarseMesh
../../_images/120c4a40121c9c865bd87a1d7e493404e6322b0103117ce78dee9338f9287ddb.png

In red, requested enriched nodes. In orange, extra enriched node related to the mesh transition

[output:0]
../../_images/05f8553eecdc4c39e0f580ed10be050b4d658f2c3c9bb20abe9b66b0ac476792.png

Show patches#

Hide code cell source

 
#nopy
topoc=cdomain.topology
adj=topoc.connectivity(0,dim)
enriched=j.getEnriched
nfh=0
nfv=0
plt=pv.Plotter(shape=(3, 3))
for i in enriched:
    k=0
    for c in adj.links(i):
        if k>0:
            fine=np.append(fine,j.getChildren(c))
        else:
            fine=j.getChildren(c)
            k=1
    if len(fine)>0 :
        tdomain =mesh.create_submesh(fdomain,dim,fine)
        tcells, ttypes, tx = plot.vtk_mesh(tdomain[0])
        tgrid = pv.UnstructuredGrid(tcells, ttypes, tx)
        if nfv>2:
            nfv=0
            plt.show()
            plt=pv.Plotter(shape=(3, 3))
            plt.subplot(0, 0)
            plt.add_mesh(tgrid,**pvopt)
            plt.add_title("Enriched node {}".format(i),font_size=7)
        else:
            plt.subplot(nfv, nfh)
            plt.add_mesh(tgrid,**pvopt)
            plt.add_title("Enriched node {}".format(i),font_size=7)
            nfh=nfh+1
            if nfh>2:
                nfv=nfv+1
                nfh=0
plt.show()
../../_images/50fd196f38f2eb482c22396469e4f75d492d8f0d7f79781f5eed6ce26594143b.png

Hide code cell source

 
#nopy
#topoc=cdomain.topology
#adj=topoc.connectivity(0,dim)
enriched=j.getExtraEnriched
nfh=0
nfv=0
plt=pv.Plotter(shape=(3, 3))
for i in enriched:
    k=0
    for c in adj.links(i):
        if k>0:
            fine=np.append(fine,j.getChildren(c))
        else:
            fine=j.getChildren(c)
            k=1
    if len(fine)>0 :
        tdomain =mesh.create_submesh(fdomain,dim,fine)
        tcells, ttypes, tx = plot.vtk_mesh(tdomain[0])
        tgrid = pv.UnstructuredGrid(tcells, ttypes, tx)
        if nfv>2:
            nfv=0
            plt.show()
            plt=pv.Plotter(shape=(3, 3))
            plt.subplot(0, 0)
            plt.add_mesh(tgrid,**pvopt)
            plt.add_title("Extra enriched node {}".format(i),font_size=7)
        else:
            plt.subplot(nfv, nfh)
            plt.add_mesh(tgrid,**pvopt)
            plt.add_title("Extra enriched node {}".format(i),font_size=7)
            nfh=nfh+1
            if nfh>2:
                nfv=nfv+1
                nfh=0
plt.show()
../../_images/692f01e65cd1c8d9d7f547b4a4d64774610bdc006a891a7a146437e5dc6e81b0.png

Show master#

[output:0]
../../_images/5e27acace07eada30044718c5d60faa8c096cc4f1747459869b2cdbc32183539.png 
[stderr:1] 2026-09-23 15:48:36.978 (   3.451s) [    7F5BDBD78140]vtkXOpenGLRenderWindow.:1450  WARN| bad X server connection. DISPLAY=

[output:1]
../../_images/7f464bae1d5c4959ef680cfce3b1d8d59ca139ddda802b550668aa35d12ab9e3.png

Damage#

damage function: hat function with damage=1 on \(x=\frac{L}{2s}(\frac{y}{B}+s)\) and damage=0 on \(x=\frac{L}{2s}(\frac{y}{B}+s ) \pm \Delta\)

%%px --local
def dam(x):
    res=np.zeros(len(x[0]))
    t1=L*(x[1]/B+s)/(2*s)
    planep=np.logical_and(np.greater(x[0],t1),np.less(x[0],t1+delta))
    res[planep]=L*x[1][planep]/(2*delta*s*B)-x[0][planep]/delta+1+L/(2*delta)
    planem=np.logical_and(np.less_equal(x[0],t1),np.greater(x[0],t1-delta))
    res[planem]=-L*x[1][planem]/(2*delta*s*B)+x[0][planem]/delta+1-L/(2*delta)
    return res

damage space and field

%%px --local
element_scal = basix.ufl.element("Lagrange", "triangle", 1)
space_scal = fem.functionspace(fdomain, element_scal)
space_scalc = fem.functionspace(cdomain, element_scal)
d=fem.Function(space_scal,name="d")
d.interpolate(dam)
dc=fem.Function(space_scalc,name="dc")
dc.interpolate(dam)
../../_images/909ab3b6f55a1fd4624aaaa6ec4f3f96d57b83f85f319ae1fbe339c5f40965d1.png

Linear problem directly on fine mesh#

Simple elasticity problem description#

Space#

%%px --local
space = fem.functionspace(fdomain, fdomain.ufl_domain().ufl_coordinate_element())

Lamé#

%%px --local
E=5.e+6
nu=0.3
mu=fem.Constant(fdomain,E / (2.0 * (1.0 + nu)))
lmbda = fem.Constant(fdomain,E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)))

Formulation#

%%px --local
u = ufl.TrialFunction(space)
v = ufl.TestFunction(space)

%%px --local
def eps(v):
    return 0.5*(ufl.grad(v) + ufl.grad(v).T)

%%px --local
def sigma(strain): 
    return (1-d)*(2.0*mu*strain + lmbda*ufl.tr(strain)*ufl.Identity(dim))

%%px --local
a_ufl=ufl.inner(sigma(eps(u)), eps(v)) * ufl.dx
f=fem.Function(space,name='load')
f.x.array[:]=0.
b_ufl=ufl.dot(f,v)*ufl.dx

Dirichlet#

%%px --local
if trac==True:
    imp=0.01
    imp2=-0.09
    imp3=0.08
else:
    imp=-0.01
    imp2=-0.09
    imp3=0.0
if D1d:
    clamp=np.array([0.])
    loading=np.array([imp])
else:
    if D2d:
        clamp=np.array([0.,0.])
        loading=np.array([imp,imp2])
    else:
        clamp=np.array([0.,0.,0.])
        loading=np.array([imp,imp2,imp3])
bc_nodes=mesh.locate_entities(fdomain, 0, lambda x: np.isclose(x[0], 0.))
bdofs=fem.locate_dofs_topological(space,0,bc_nodes)
bcb=fem.dirichletbc(clamp,bdofs,space)
imp_nodes=mesh.locate_entities(fdomain, 0, lambda x: np.isclose(x[0], L))
idofs=fem.locate_dofs_topological(space,0,imp_nodes)
bci=fem.dirichletbc(loading,idofs,space)
bcs=[bcb,bci]

Field#

%%px --local
disp=fem.Function(space,name='disp_fine')

Generate Multipoint constraint#

%%px --local
mpc=core.generateMPC(j,space)

Solve without mpc#

%%px --local
problem = petsc.LinearProblem(a_ufl, b_ufl,bcs=bcs,petsc_options=petsc_options,petsc_options_prefix="lpnmpc_")

%%px --local
disp = problem.solve()

Solve with mpc#

%%px --local
problemw = LinearProblem_mpc(a_ufl, b_ufl,mpc,bcs=bcs,petsc_options=petsc_options,petsc_options_prefix="lpmpc_")

%%px --local
dispw = problemw.solve()

plot displacement solution at fine scale with or without MPC#

../../_images/c39a390a8ccd29b8bdf9e07f7afd357d4cb3b4a3248f4e275f5613a84cb94725.png 
[output:0]
../../_images/96a2fc19cde7341ead6cf66323526d5a118bc39ceff7eb11259f8dbe400dd899.png

Two-scale approach#

Use nested strategy

Enriched space#

%%px --local

el_mixed = basix.ufl.mixed_element([cdomain.ufl_domain().ufl_coordinate_element(), cdomain.ufl_domain().ufl_coordinate_element()])
coarse_enriched_space = fem.functionspace(cdomain, el_mixed)

std_=coarse_enriched_space.sub(0)
std, std_to_mix=std_.collapse()
enr_=coarse_enriched_space.sub(1)
enr, enr_to_mix=enr_.collapse()
std_to_mix_np=np.array(std_to_mix[0],dtype=int)

disp_ce=fem.Function(coarse_enriched_space,name='disp_coarse_enriched')

Solve coarse non enriched problem (standard part) to init two scale loop#

Formulation#

%%px --local

muc=fem.Constant(cdomain,E / (2.0 * (1.0 + nu)))
lmbdac = fem.Constant(cdomain,E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)))
ucs = ufl.TrialFunction(std)
vcs = ufl.TestFunction(std)
def sigmac(strain): 
    return (1-dc)*(2.0*muc*strain + lmbdac*ufl.tr(strain)*ufl.Identity(dim))
acs=ufl.inner(sigmac(eps(ucs)), eps(vcs)) * ufl.dx
fcs=fem.Function(std,name='load')
fcs.x.array[:]=0.
bcst=ufl.dot(fcs,vcs)*ufl.dx

Boundary condition#

%%px
bc_nodesc=mesh.locate_entities(cdomain, 0, lambda x: np.isclose(x[0], 0.))
bdofscs=fem.locate_dofs_topological(std,0,bc_nodesc)
bcbcs=fem.dirichletbc(clamp,bdofscs,std)
imp_nodesc=mesh.locate_entities(cdomain, 0, lambda x: np.isclose(x[0], L))
idofscs=fem.locate_dofs_topological(std,0,imp_nodesc)
bcics=fem.dirichletbc(loading,idofscs,std)
bcscs=[bcbcs,bcics]

Resolution#

%%px
problemcs = petsc.LinearProblem(acs, bcst,bcs=bcscs,petsc_options=petsc_options,petsc_options_prefix="lptsc_")
dispcs = problemcs.solve()

[output:0]
../../_images/0ed9867a2fb2f165f80c6ec7ae706926432c911bb602eb9788980b3c0568054e.png

Store coarse standard in coarse enriched field#

%%px 
if nested_strategy:
    disp_ce=dispcs
else:
    disp_ce.x.array[std_to_mix_np]=dispcs.x.array
#print(disp_ce.x.array)

Twoscale loop#

Choose an enriched function#

%%px --local
enriched_shift=core.generateEnrichedShiftFunction(disp)

Resolution with linearBasicLoop (i.e. scale loop resolution)#

Dirichlet boundary condition must be defined on mixed space and merged into one Dirichlet BC object

%%px 
if nested_strategy:
    BC_ce=bcscs
else:
    bc_nodesc=mesh.locate_entities(cdomain, 0, lambda x: np.isclose(x[0], 0.))
    imp_nodesc=mesh.locate_entities(cdomain, 0, lambda x: np.isclose(x[0], L))
    clamp_std_dofs=fem.locate_dofs_topological(std_,0,bc_nodesc)
    all_BC_values=fem.Function(coarse_enriched_space)
    imp_std_dof_x=fem.locate_dofs_topological(std.sub(0),0,imp_nodesc)
    imp_std_dof_x_ce=std_to_mix_np[imp_std_dof_x]
    if D1d:
       all_BC_dof=np.concat((clamp_std_dofs,imp_std_dof_x_ce))
       all_BC_values.x.array[imp_std_dof_x_ce]=imp
    else:
       if D2d:
          imp_std_dof_y=fem.locate_dofs_topological(std.sub(1),0,imp_nodesc)
          imp_std_dof_y_ce=std_to_mix_np[imp_std_dof_y]
          all_BC_dof=np.concat((clamp_std_dofs,imp_std_dof_x_ce,imp_std_dof_y_ce))
          all_BC_values.x.array[imp_std_dof_x_ce]=imp
          all_BC_values.x.array[imp_std_dof_y_ce]=imp2
       else:
          imp_std_dof_y=fem.locate_dofs_topological(std.sub(1),0,imp_nodesc)
          imp_std_dof_y_ce=std_to_mix_np[imp_std_dof_y]
          imp_std_dof_z=fem.locate_dofs_topological(std.sub(2),0,imp_nodesc)
          imp_std_dof_z_ce=std_to_mix_np[imp_std_dof_z]
          all_BC_dof=np.concat((clamp_std_dofs,imp_std_dof_x_ce,imp_std_dof_y_ce,imp_std_dof_z_ce))
          all_BC_values.x.array[imp_std_dof_x_ce]=imp
          all_BC_values.x.array[imp_std_dof_y_ce]=imp2
          all_BC_values.x.array[imp_std_dof_z_ce]=imp3
    all_BC_dof.sort()
    BC_ce=[fem.dirichletbc(all_BC_values,all_BC_dof)]

Systems are created and TS resolution is done:

%%px 
#--local
[A,AD,BND,BD]=util.createFineScaleSytems(a_ufl,b_ufl,bcs=bcs,MPC=mpc)
[dispf, r,nm, it,hrb,hrr]=linear.linearBasicLoop(j,dispw.function_space,A,AD,BND,BD,enriched_shift,disp_ce,mpc,BC_ce,itmax,epsr)

Solutions with direct solver and TS solver are plotted bellow with their differences

[output:0]
../../_images/aed0a4ccb08522915ba27cc71f7d2d4f085cd2932b75a8a60df7604396df8e48.png 
[output:0]
../../_images/8df42286a011847a6d9f26b277d2ac7131ddfa8ca94aec101e9350d5b67caf51.png

The following residual evolution histories are plotted below:

  • TS residual fine system relative to rhs \(\frac{|AD.S_{ts_i}-BD|}{|BD|}\)

  • TS residual fine system relative to first residual \(\frac{|AD.S_{ts_i}-BD|}{|AD.S_{ts_0}-BD|}\)

  • TS relative solution difference evolution \(\frac{|S_{ts_i}-S_{ts_{i-1}})|}{|S_{ts_0}|}\)

Text(0, 0.5, 'Criterion')
../../_images/e11989cc5c9b5cf9f9e1ecaf8346387d47ee3e476d94b7b800696dfe4bd4c2d6.png

The TS solver solution iterations are:

../../_images/local_damage.gif 
Text(0, 0.5, 'Criterion')
../../_images/b363e46a644d85654fee3bd66bce736024ccdf614b92919d40e3335699ebb70f.png