Pure discontinuous test case

Pure discontinuous test case#

Initialize#

For parallel testing#

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

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:55.583 (   1.028s) [    7F397B9F7140]vtkXOpenGLRenderWindow.:1450  WARN| bad X server connection. DISPLAY=
../../_images/aaf1301aa8559729a22a77a89e19beae2a75628b6350b967e96733765c98fc86.png 
[stderr:0] 2026-09-23 15:48:55.826 (   1.409s) [    7FFB942A1140]vtkXOpenGLRenderWindow.:1450  WARN| bad X server connection. DISPLAY=

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

Choose a mesh refinement localization criterion#

Finest mesh close to the middle. The middle is fixed by posmx and is shifted from the true middle position \(\frac{L}{2}\)

%%px --local
posmx=L/2+L/(10*s)

def crit1(coords,level):
    return np.logical_and(coords[0]<posmx+L/(8*s),coords[0]>posmx-L/(8*s))

Choose macro node to enrich#

select central region around \(x=\frac{L}{2} \pm \frac{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/e18efae6d56becd4f2353f720611c847cc13e9454fae8a3a5741203508ffba33.png

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

[output:0]
../../_images/83af26dbf0468b3d9ed565ea3894dc00674be6cdcd59fcdc5a3883774705b808.png

Damage#

A damage field is used to impose a pure discontinuous displacement in the middle (i.e. posmx). In between \(x=posmx \pm \Delta\), damage is one and is null otherwise. Thus only the fine scale is able to get an element fully damaged and have discontinuous behavior.

%%px --local
delta=2*L/(7*s*ref)
def dam(x):
    res=np.zeros(len(x[0]))
    disco=np.logical_and(np.greater(x[0],posmx-delta),np.less(x[0],posmx+delta))
    res[disco]=0.999
    return res

damage space and field

%%px --local
if D2d:
    element_scal = basix.ufl.element("Lagrange", "triangle", 1)
else:
   element_scal = basix.ufl.element("Lagrange", "tetrahedron", 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)

[output:0]
../../_images/8c023d6e7ae8ce15b687c467fed477d4c59ed639504f05956f25cf12427a9e1e.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 D1d:
    clamp=np.array([0.])
    loading=np.array([imp])
else:
    if D2d:
        clamp=np.array([0.,0.])
        loading=np.array([imp,0.])
    else:
        clamp=np.array([0.,0.,0.])
        loading=np.array([imp,0.,0.])
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="f_")

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

Solve with mpc#

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

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

Plot displacement solution at fine scale with or without MPC#

../../_images/d03a207af7d6553250436bfe666eb591dd00e450d2d8e9c4197b95848fa25308.png 
[output:0]
../../_images/a6b9ac71c529c1b7dc70375a02ce9197166e9d421470f8c29ced03c09e0b8f01.png

Two-scale approach#

[stdout:0] 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 initialize 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="c_")
dispcs = problemcs.solve()

[output:0]
../../_images/70c5ffed30a3ffd570750c55a1133c3e32bae507e5d93cfb37c05ef9ae151fad.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)

Two-scale 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
        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_dof.sort()
    BC_ce=[fem.dirichletbc(all_BC_values,all_BC_dof)]

Systems are created and TS resolution is done:

%%px 
[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 below with their differences

[output:0]
../../_images/2230e193f291692aa4a8378d8ee9c724d18cf936433db5ed1277276bb52a3793.png 
[output:0]
../../_images/d23d7f09a3dc32d6699b201ed66c42d748488218dd8b5b4d070c6f836ff83eed.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}|}\)

[stdout:0] At twoscale iteration 0 residual is 0.06601264096891646(b) 1.(r) 1.(ds)
At twoscale iteration 1 residual is 0.22761885671455576(b) 3.448110140325021(r) 0.6179566252624156(ds)
At twoscale iteration 2 residual is 0.040770278648964965(b) 0.6176132033281713(r) 0.18979175531310466(ds)
At twoscale iteration 3 residual is 0.0036263635198145227(b) 0.05493438024274893(r) 0.02541356553501255(ds)
At twoscale iteration 4 residual is 0.0008415769628363971(b) 0.012748724342549372(r) 0.003911985644169589(ds)
At twoscale iteration 5 residual is 0.0005623952684586594(b) 0.008519508691122902(r) 0.0017640981704001388(ds)
At twoscale iteration 6 residual is 0.000458440379444661(b) 0.0069447362310580485(r) 0.001240441362790908(ds)
At twoscale iteration 7 residual is 0.0003831495285991963(b) 0.005804184213438922(r) 0.0009500587764458796(ds)
At twoscale iteration 8 residual is 0.00032139987424050143(b) 0.004868762550976256(r) 0.000756366691911035(ds)
At twoscale iteration 9 residual is 0.0002715758355946261(b) 0.004113997434559598(r) 0.0006148419260741672(ds)
At twoscale iteration 10 residual is 0.00022770890575930238(b) 0.0034494742585215498(r) 0.0005144738250480318(ds)
At twoscale iteration 11 residual is 0.00019220530262613225(b) 0.0029116438882764357(r) 0.00042498863128257033(ds)
At twoscale iteration 12 residual is 0.00016898581638202854(b) 0.0025599008599216525(r) 0.0003521594676843733(ds)
At twoscale iteration 13 residual is 0.0001414288849649503(b) 0.0021424515500227496(r) 0.00030533165413107097(ds)
At twoscale iteration 14 residual is 0.00012187750801286742(b) 0.0018462752924891493(r) 0.0002552269363930402(ds)
At twoscale iteration 15 residual is 0.00011037675615855308(b) 0.0016720548449277535(r) 0.0002059197667283898(ds)
At twoscale iteration 16 residual is 9.178574317545515e-05(b) 0.0013904267702101866(r) 0.00018920059977387278(ds)
At twoscale iteration 17 residual is 7.756263178517711e-05(b) 0.0011749663495768824(r) 0.0001630947876054393(ds)
At twoscale iteration 18 residual is 6.223713559078065e-05(b) 0.0009428063273530657(r) 0.00014494403229445342(ds)
At twoscale iteration 19 residual is 5.3837121722810284e-05(b) 0.000815557761855956(r) 0.00011639341359681859(ds)
At twoscale iteration 20 residual is 4.721110933634721e-05(b) 0.0007151828595765111(r) 0.0001005300685767685(ds)
At twoscale iteration 21 residual is 4.6512587114274845e-05(b) 0.0007046012162454817(r) 8.707613634518074e-05(ds)
At twoscale iteration 22 residual is 3.4211803812476115e-05(b) 0.0005182614013062365(r) 8.190682452019033e-05(ds)
At twoscale iteration 23 residual is 2.852407415761205e-05(b) 0.00043210018170676214(r) 6.44481613494334e-05(ds)
At twoscale iteration 24 residual is 2.431147292195804e-05(b) 0.0003682851127469002(r) 5.4923158210260494e-05(ds)
At twoscale iteration 25 residual is 2.113895950586614e-05(b) 0.00032022593242133577(r) 4.114620036313281e-05(ds)
At twoscale iteration 26 residual is 2.2332927472664827e-05(b) 0.0003383128919683851(r) 4.6868558719009934e-05(ds)
At twoscale iteration 27 residual is 1.961986924365145e-05(b) 0.00029721382080274457(r) 3.461357060099731e-05(ds)
At twoscale iteration 28 residual is 1.40370796442118e-05(b) 0.0002126422975687562(r) 2.8080280483168364e-05(ds)
At twoscale iteration 29 residual is 2.3987022838403934e-05(b) 0.0003633701437532058(r) 6.202493313492815e-05(ds)
At twoscale iteration 30 residual is 1.6564663091269924e-05(b) 0.000250931682904033(r) 3.536332346244362e-05(ds)

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

The TS solver solution iterations are:

../../_images/pure_discontinuous.gif 
[stdout:0] With the nested strategy, convergence is not reached either with TS used as solver or as preconditioner.
It is certainly related to the iterative solver at the coarse level, perturbed by poor conditioning.
With monolithic approach, coarse solver is a LU direct solver and thus is less perturbed.
Clearly, this test case would be better if no coarse enriched nodes exist, and only extra enrichment nodes exist.
In this case, no patch would be close to rigid body mode movement...