Analytical 3D elasticity

Analytical 3D elasticity#

Setting#

For parallel testing#

Create a server controlling n>1 mpi engines (consider openmpi as mpi backend)

Starting 2 engines with <class 'ipyparallel.cluster.launcher.MPIEngineSetLauncher'>

  0%|          | 0/2 [00:00<?, ?engine/s]

 50%|█████     | 1/2 [00:05<00:05,  5.68s/engine]

100%|██████████| 2/2 [00:05<00:00,  2.84s/engine]

[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 twosacle library#

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

Configuration#

Basic test case#

[stderr:0] 2026-09-23 15:49:19.064 (   1.137s) [    7FD766E0F140]vtkXOpenGLRenderWindow.:1450  WARN| bad X server connection. DISPLAY=

[output:0]

../../_images/97fb60e0015de4fb82a22ab2e50de5977eca80e48eda0563e24bc8009fb8890b.png

Choose macro node to enrich#

Select all nodes

%%px
def enriched(coords):
    return coords[0]>-1.

Scale Jump construction#

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

%%px
fdomain=j.getFineMesh
cdomain=j.getCoarseMesh

[output:0]

../../_images/46d62e47e7d5f04afb3626ef07df9cf9b2ce8411b33f3ae2444e2393393735ce.png

Fine-scale discretisation#

Space#

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

Lamé#

%%px
E=36.5e+9
nu=0.2
F=1.e+10
mu=fem.Constant(fdomain,E / (2.0 * (1.0 + nu)))
lmbda = fem.Constant(fdomain,E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)))
coef1=(nu+1)*(1-2.*nu)*F/E
coef2=(1.-nu)*F

Exact solution#

%%px
def exact_sol(coord):
    x=coef1*(2.*nu*(np.power(coord[1],2)-np.power(coord[2],2))+np.power(coord[0],2)*(L/2-coord[0]/3.))
    y=-4*coef1*nu*coord[0]*coord[1]
    z=4*coef1*nu*coord[0]*coord[2]
    return np.array([x,y,z])

%%px
disp_exact=fem.Function(space,name='disp_exact')
disp_exact.interpolate(exact_sol)

[output:0]

../../_images/d886c5a53a6f161a59a59995e3c14ab35ca666bb70293a1d48e928165df61f41.png

Formulation#

Bilinear form#

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

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

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

%%px
a_ufl=ufl.inner(sigma(eps(u)), eps(v)) * ufl.dx

Volume load#

%%px
metadatav = {"quadrature_rule": "default", "quadrature_degree": 2}
dxx = ufl.measure.Measure("dx", domain = fdomain, metadata = metadatav)
x,y,z = ufl.SpatialCoordinate(fdomain)
fv=ufl.as_vector([coef2*(2*x-L),0,0])
b_ufl=ufl.dot(fv,v)*dxx

Surface load#

Define each face to integrate for surface loading.

%%px
bc_face_xz0=mesh.locate_entities(fdomain, 2, lambda x: np.isclose(x[1], 0.))
bc_face_xzB=mesh.locate_entities(fdomain, 2, lambda x: np.isclose(x[1], B))
bc_face_xy0=mesh.locate_entities(fdomain, 2, lambda x: np.isclose(x[2], 0.))
bc_face_xyH=mesh.locate_entities(fdomain, 2, lambda x: np.isclose(x[2], H))
facet_indices=np.hstack((bc_face_xz0,bc_face_xzB,bc_face_xy0,bc_face_xyH)).astype(np.int32)
facet_markers = np.hstack((np.full_like(bc_face_xz0, 1),np.full_like(bc_face_xzB, 1),np.full_like(bc_face_xy0, 2),np.full_like(bc_face_xyH, 2))).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = mesh.meshtags(fdomain, dim, facet_indices[sorted_facets], facet_markers[sorted_facets])
metadatas = {"quadrature_rule": "default", "quadrature_degree": 3}
ds = ufl.Measure("ds", domain=fdomain, subdomain_data=facet_tag,metadata=metadatas) 

Define each formulation based on a UFL expression for each surface.

%%px
n=ufl.FacetNormal(fdomain)
fsxy=x*F*nu*(L-x+4*(1-2*nu))
fsxz=x*F*nu*(L-x-4*(1-2*nu))
b_ufl+=ufl.inner(fsxz*n,v)*ds(1)
b_ufl+=ufl.inner(fsxy*n,v)*ds(2)
#print(b_ufl)

Dirichlet#

A(0,0,0) corner clamped to block x,y,z translations and \(\theta_x\),\(\theta_y\),\(\theta_z\) rotations:

%%px
node_A=mesh.locate_entities(fdomain, 0, lambda x: np.logical_and(np.isclose(x[0], 0.),np.logical_and(np.isclose(x[1], 0.),np.isclose(x[2],0.))))
A_dofs=fem.locate_dofs_topological(space,0,node_A)
bcA=fem.dirichletbc(np.array([0.,0.,0.]),A_dofs,space)

B(L,0,0) node clamped in the y,z directions to block \(\theta_z\) and \(\theta_y\):

%%px
node_B=mesh.locate_entities(fdomain, 0, lambda x: np.logical_and(np.isclose(x[0], L),np.logical_and(np.isclose(x[1], 0.),np.isclose(x[2],0.))))
By_dofs=fem.locate_dofs_topological(space.sub(1),0,node_B)
bcBy=fem.dirichletbc(0.,By_dofs,space.sub(1))
Bz_dofs=fem.locate_dofs_topological(space.sub(2),0,node_B)
bcBz=fem.dirichletbc(0.,Bz_dofs,space.sub(2))

C(0,B,0) node clamped in the z direction to block \(\theta_x\)

%%px
node_C=mesh.locate_entities(fdomain, 0, lambda x: np.logical_and(np.isclose(x[0], 0.),np.logical_and(np.isclose(x[1], B),np.isclose(x[2],0.))))
C_dofs=fem.locate_dofs_topological(space.sub(2),0,node_C)
bcC=fem.dirichletbc(0.,C_dofs,space.sub(2))

bcs=[bcA,bcBy,bcBz,bcC]

Field#

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

Solve with direct resolution#

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

%%px
disp = problem.solve()

[output:0]

../../_images/314ed1b01aad7d15780cfa32e34f0284883f41e00ca03e4dfd77c019036b215f.png

Two-scale approach#

[stdout:0] Use nested strategy

Enriched space#

%%px
if nested_strategy:
    std=fem.functionspace(cdomain, cdomain.ufl_domain().ufl_coordinate_element())
else:
    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

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 2.0*muc*strain + lmbdac*ufl.tr(strain)*ufl.Identity(dim)
acs=ufl.inner(sigmac(eps(ucs)), eps(vcs)) * ufl.dx
dxxc = ufl.measure.Measure("dx", domain = cdomain, metadata = metadatav)
xc,_,_ = ufl.SpatialCoordinate(cdomain)
fvc=ufl.as_vector([coef2*(2*xc-L),0,0])
bcst=ufl.dot(fvc,vcs)*dxxc
bc_face_xz0c=mesh.locate_entities(cdomain, 2, lambda x: np.isclose(x[1], 0.))
bc_face_xzBc=mesh.locate_entities(cdomain, 2, lambda x: np.isclose(x[1], B))
bc_face_xy0c=mesh.locate_entities(cdomain, 2, lambda x: np.isclose(x[2], 0.))
bc_face_xyHc=mesh.locate_entities(cdomain, 2, lambda x: np.isclose(x[2], H))
facet_indicesc=np.hstack((bc_face_xz0c,bc_face_xzBc,bc_face_xy0c,bc_face_xyHc)).astype(np.int32)
facet_markersc = np.hstack((np.full_like(bc_face_xz0c, 1),np.full_like(bc_face_xzBc, 1),np.full_like(bc_face_xy0c, 2),np.full_like(bc_face_xyHc, 2))).astype(np.int32)
sorted_facetsc = np.argsort(facet_indicesc)
facet_tagc = mesh.meshtags(cdomain, dim, facet_indicesc[sorted_facetsc], facet_markersc[sorted_facetsc])
dsc = ufl.Measure("ds", domain=cdomain, subdomain_data=facet_tagc,metadata=metadatas) 
nc=ufl.FacetNormal(cdomain)
fsxyc=xc*F*nu*(L-xc+4*(1-2*nu))
fsxzc=xc*F*nu*(L-xc-4*(1-2*nu))
bcst+=ufl.inner(fsxzc*nc,vcs)*dsc(1)
bcst+=ufl.inner(fsxyc*nc,vcs)*dsc(2)

Boundary condition#

%%px
node_Ac=mesh.locate_entities(cdomain, 0, lambda x: np.logical_and(np.isclose(x[0], 0.),np.logical_and(np.isclose(x[1], 0.),np.isclose(x[2],0.))))
A_dofsc=fem.locate_dofs_topological(std,0,node_Ac)
bcAc=fem.dirichletbc(np.array([0.,0.,0.]),A_dofsc,std)
node_Bc=mesh.locate_entities(cdomain, 0, lambda x: np.logical_and(np.isclose(x[0], L),np.logical_and(np.isclose(x[1], 0.),np.isclose(x[2],0.))))
By_dofsc=fem.locate_dofs_topological(std.sub(1),0,node_Bc)
bcByc=fem.dirichletbc(0.,By_dofsc,std.sub(1))
Bz_dofsc=fem.locate_dofs_topological(std.sub(2),0,node_Bc)
bcBzc=fem.dirichletbc(0.,Bz_dofsc,std.sub(2))
node_Cc=mesh.locate_entities(cdomain, 0, lambda x: np.logical_and(np.isclose(x[0], 0.),np.logical_and(np.isclose(x[1], B),np.isclose(x[2],0.))))
C_dofsc=fem.locate_dofs_topological(std.sub(2),0,node_Cc)
bcCc=fem.dirichletbc(0.,C_dofsc,std.sub(2))

bcscs=[bcAc,bcByc,bcBzc,bcCc]

Resolution#

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

[output:0]

../../_images/c4943d993f58510d116bd0f09fa50c6402dfa29a9e5864705a154647272ec3d8.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
enriched_shift=core.generateEnrichedShiftFunction(disp)

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

Dirichlet boundary conditions must be defined on the mixed space and merged into one Dirichlet BC object for non-nested strategy.

%%px 
if nested_strategy:
    BC_ce=bcscs
else:
    A_dofsce=fem.locate_dofs_topological(std_,0,node_Ac)
    By_dofsce=fem.locate_dofs_topological(std.sub(1),0,node_Bc)
    By_dofsce=std_to_mix_np[By_dofsce]
    Bz_dofsce=fem.locate_dofs_topological(std.sub(2),0,node_Bc)
    Bz_dofsce=std_to_mix_np[Bz_dofsce]
    C_dofsce=fem.locate_dofs_topological(std.sub(2),0,node_Cc)
    C_dofsce=std_to_mix_np[C_dofsce]
    all_BC_dof=np.concat((A_dofsce,By_dofsce,Bz_dofsce,C_dofsce))
    all_BC_dof.sort()
    all_BC_values=fem.Function(coarse_enriched_space)
    BC_ce=[fem.dirichletbc(all_BC_values,all_BC_dof)]

Systems are created and TS solution is computed considering Dirichlet BC (i.e. rigide body mode elimination) imposed at patch level:

%%px 
[A,AD,BND,BD]=util.createFineScaleSytems(a_ufl,b_ufl,bcs=bcs)
[dispf, r, nm,it,hr,hrr]=linear.linearBasicLoop(j,space,A,AD,BND,BD,enriched_shift,disp_ce,None,BC_ce,itmax*2,epsr)

Resolution with init/runLinearBasicLoop (i.e. scale loop resolution)#

When Dirichlet BC are only used to block rigide body mode, patch problems do not need to block these DOFs. Because at patch level blocking these DOFs is meaningless and over constrains the patch. To avoid using these DOFs one must use the init/runLinearBasicLoop function where system provided to initLinearBasicLoop is the one without Dirichlet BC eliminated.

%%px 
[dispf2, cm, pm]=linear.initLinearBasicLoop(j,A,BND,space,disp_ce,None,BC_ce)
[r2, it2,hr2,hrr2]=linear.runLinearBasicLoop(dispf2,cm,pm,A,AD,BND,BD,enriched_shift,disp_ce,None,itmax,epsr)

Convergence to 1e-07 obtained after  20 TS iterations with relative convergence 6.409758058093389e-08

resolution with init/runLinearBasicLoop with Eimp#

%%px 
if nested_strategy:
    [dispf3, cm3, pm3]=linear.initLinearBasicLoop(j,A,BND,space,disp_ce,None,BC_ce)
    [r3, it3,hr3,hrr3]=linear.runLinearBasicLoopEImp(dispf3,cm3,pm3,A,AD,BND,BD,enriched_shift,disp_ce,None,5,itmax,epsr)

Convergence to 1e-07 obtained after  22 TS iterations with relative convergence 7.37988951083958e-08

Results#

Solutions with the direct solver and TS solver are plotted bellow, together with their differences

[output:0]

../../_images/7fe14a19d50527c8a052897b283e7b9db0e3928628057d77751bc93d20887387.png

[output:0]

../../_images/307a629265b2f7f9a04f14021877a7c168baab0f2226446382935ddb3fbc0872.png

[output:0]

../../_images/97a32124afa812b35c80dfafea0aa2a782ba3949cac3ef3aa9678fc2272e3271.png

[output:0]

../../_images/7289259c8612353bce5eb4643fdb528da871f182b3143b18f830f624fdb11772.png

[output:0]

../../_images/029bfed3e4ba6ab12953347b81c8a451bdbcfdc720b267665b8e0a3d5ab6628d.png

[output:0]

../../_images/dcc0ecfcf22f93b1338c3a499bf8e96376390e179858431a48ee95d043c3ddc5.png

The difference between imposing and not Dirichlet BC at patch level may looks almost “null” as it can be seen above. But depending on scale jump TS can not converge with Dirichlet imposed at patch level. In all cases, the way to reach this solution is different as it can be seen below with residual evolution:

Text(0, 0.5, '$norm(AD.S_{ts_i}-BD)/norm(BD)$')

../../_images/700308a478fa734b9f47b4d47e6cb10ce3aa1dbfa355235fd826ba28aa185525.png

The following criteria have been also computed manually:

relative difference TS-direct \(\frac{|AD.(S_{ts_i}-S_{direct})|}{|AD.S_{exact}|}\)
relative difference TS-exact \(\frac{|AD.(S_{ts_i}-S_{exact})|}{|AD.S_{exact}|}\)
relative difference direct-exact \(\frac{|AD.(S_{direct}-S_{exact})|}{|AD.S_{exact}|}=cst\)

Text(0, 0.5, 'Criterion')

../../_images/30a5ed2d1a5258ac7311f017dbb244c5493972968af139e9e92cdecef789f93f.png

Analytical 3D elasticity

Contents

Analytical 3D elasticity#

Setting#

For parallel testing#

Import twosacle library#

Configuration#

Basic test case#

Choose a mesh refinement localization criterion#

Choose macro node to enrich#

Scale Jump construction#

Fine-scale discretisation#

Space#

Lamé#

Exact solution#

Formulation#

Bilinear form#

Volume load#

Surface load#

Dirichlet#

Field#

Solve with direct resolution#

Two-scale approach#

Enriched space#

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

Formulation#

Boundary condition#

Resolution#

Store coarse standard in coarse enriched field#

Twoscale loop#

Choose an enriched function#

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

Resolution with init/runLinearBasicLoop (i.e. scale loop resolution)#

resolution with init/runLinearBasicLoop with Eimp#

Results#