Adapt MeshFem with XFEM enrichment

Questions and issues
1

Hello,

I am working on a model that uses level sets to define an interface that moves through time with some enrichment around that interface. It is straightforward to implement this for a single time step (e.g. the demo_crack example). However, I need to be able to update the location of the interface.

In the case where the only enrichment is the heaviside function using mesh_fem_level_set, e.g.:

mf_u = gf.MeshFem(‘levelset’,mls,mf_pre_u)

this can easily be adapted along with the mesh and integration methods.

However, if there is an additional enrichment, such as the use of the cutoff function, in which the MeshFem is the product of the heaviside-enriched MeshFem and an additional global function, e.g.:

mfls_u = gf.MeshFem(‘levelset’,mls,mf_pre_u)
mf_sing_u = gf.MeshFem(‘global function’,m,ls,[ckoff0,ckoff1,ckoff2,ckoff3],1)
mf_u = gf.MeshFem(‘sum’,mf_sing_u,mfls_u)

it does not seem possible to adapt the final MeshFem through the Python interface. As a result, the second solve will fail. Is there something I am missing here or a simple way around this issue?
I have a simple reproducible example below based on the demo_crack example (here I simply update the values using another level set but of course in my full model I update based on an evolution law).

Thank you for your kind attention,
Janine Birnbaum

Parameters:

nx = 20 #
DIRICHLET = 101 #
Lambda = 1.25E10 # Lame coefficient #
Mu = 1.875E10 # Lame coefficient #
#####################################

Mesh definition:

m = gf.Mesh(‘regular_simplices’, -0.5+np.arange(nx+1)/float(nx),
-0.5+np.arange(nx+1)/float(nx))

Boundary set:

m.set_region(DIRICHLET, m.outer_faces())

Global functions for asymptotic enrichment:

ck0 = gf.GlobalFunction(‘crack’,0)
ck1 = gf.GlobalFunction(‘crack’,1)
ck2 = gf.GlobalFunction(‘crack’,2)
ck3 = gf.GlobalFunction(‘crack’,3)

coff = gf.GlobalFunction(‘cutoff’,2,0.4,0.01,0.4)
ckoff0 = ck0coff # gf.GlobalFunction(‘product’,ck0,coff)
ckoff1 = ck1coff
ckoff2 = ck2coff
ckoff3 = ck3coff

Levelset definition:

ls = gf.LevelSet(m,1,‘y’,‘x’)

mls = gf.MeshLevelSet(m)
mls.add(ls)
mls.adapt()

Basic mesh_fem without enrichment:

mf_pre_u = gf.MeshFem(m)
mf_pre_u.set_fem(gf.Fem(‘FEM_PK(2,1)’))

Definition of the enriched finite element method (MeshFemLevelSet):

mfls_u = gf.MeshFem(‘levelset’,mls,mf_pre_u)

MeshFemGlobalFunction:

mf_sing_u = gf.MeshFem(‘global function’,m,ls,[ckoff0,ckoff1,ckoff2,ckoff3],1)
mf_u = gf.MeshFem(‘sum’,mf_sing_u,mfls_u)
mf_u.set_qdim(2)

MeshIm definition (MeshImLevelSet):

mim = gf.MeshIm(‘levelset’, mls, ‘all’,
gf.Integ(‘IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(6),3)’),
gf.Integ(‘IM_STRUCTURED_COMPOSITE(IM_GAUSS_PARALLELEPIPED(2,6),9)’),
gf.Integ(‘IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(6),5)’))

Exact solution for a single crack:

mf_ue = gf.MeshFem(‘global function’,m,ls,[ck0,ck1,ck2,ck3])
A = 2+2Mu/(Lambda+2Mu); B=-2*(Lambda+Mu)/(Lambda+2*Mu)
Ue = np.zeros([2,4])
Ue[0,0] = 0; Ue[1,0] = A-B # sin(theta/2)
Ue[0,1] = A+B; Ue[1,1] = 0 # cos(theta/2)
Ue[0,2] = -B; Ue[1,2] = 0 # sin(theta/2)*sin(theta)
Ue[0,3] = 0; Ue[1,3] = B # cos(theta/2)cos(theta)
Ue /= 2np.pi
Ue = Ue.T.reshape(1,8)

Model definition:

md = gf.Model(‘real’)
md.add_fem_variable(‘u’, mf_u)

data

md.add_initialized_data(‘lambda’, [Lambda])
md.add_initialized_data(‘mu’, [Mu])
md.add_isotropic_linearized_elasticity_brick(mim,‘u’,‘lambda’,‘mu’)
md.add_initialized_fem_data(“DirichletData”, mf_ue, Ue)
md.add_Dirichlet_condition_with_penalization(mim,‘u’, 1e12, DIRICHLET, ‘DirichletData’)

Assembly of the linear system and solve:

md.solve()
U = md.variable(‘u’)

Modify the location of the level set and adapt the mesh

ls2 = gf.LevelSet(m,1,‘y’,‘x-0.01’)
ls.set_values(ls2.values(0),ls2.values(1))

mls.adapt()
mf_u.adapt()
mim.adapt()

md.solve()

1 Like
2

As far as I could check, it seems you are right, the Python interface does not support to use the adapt() method on a mesh_fem_sum that depends on mesh_fem_levelset objects, actually not even for a mesh_fem_levelset object created with
mfls2 = gf.MeshFem(“clone”, mfls1)

But this seems also to be easy to fix, I will commit some changes to the getfem git repository for fixing this limitation.

The next question will be what will happen to the variables stored inside an existing model that depend on the adapted mesh_fem. Since the number and numbering of degrees of freedom will change, any previous assigned values to the respective model variables will be lost.

For transferring fields between the non-adapted and the adapted version of a mesh_fem, one should manually do the following:

  • Create a clone copy (call it mfcopy) of the mesh_fem to be adapted (call it mf) and store the vector corresponding to the variable of interest described by this mesh_fem (call it Vcopy).
  • Adapt the mesh_fem used in the model with mf.adapt()
  • Interpolate between mfcopy and mf: Vcopy → V, using the compute_interpolate_on function ( V=gf. compute_interpolate_on(mfcopy, Vcopy, mf) )
  • Assign the computed vector to the affected variable in the model using md.set_variable(“…”, V)

AFAIK this procedure needs currently to be done by the user following these steps manually. Of course the best solution would be if we implemented an adapt() function to the getfem model, that does all this for you automatically.

1 Like
3

After some changes in GetFEM which I still have to upload to the git repository, I can run the following example (which still doesn’t use mesh_fem_sum). I need to do some more testing, also with mesh_fem_sum, and some cleaning up and upload the changes.

growing_crack

import sys
import numpy as np
import getfem as gf

np.set_printoptions(threshold=sys.maxsize)
gf.util_trace_level(1)

# Parameters:  #########################################################
E0 = 5E2            # [N/mm^2] Young modulus at reference temperature
nu = 0.45           # [-] Poisson ratio
beta = 8e-5         # [1/K] Expansion coefficient
gamma = 0.1         # [1/K] Coefficient of thermal weakening
theta0 = 30.        # [C] Reference temperature
alpha = 0.1         # [mm^2/s] Thermal diffusivity (0.1 for rubber)
alpha_env = 0.02    # [N/mm/s/K] Heat transfer coefficient
alpha_crack = 0.01  # [N/mm/s/K] Heat transfer coefficient
mR = 287e3*5e-9     # [N*mm/K] --> 287000 [N*mm/kg/K] * mass [kg]
room_temp = 30.     # [C] temperature on one side of the block
jump_temp = 100.    # [C] temperature jump between the two sides

# Mesh definition:
L = 40.             # [mm]
l = 10.             # [mm]
ny = 10             # Number of elements in each direction
quad = True         # quad mesh or triangle one

h = l/ny;
nx = np.floor(L/h);

if (quad):
  m=gf.Mesh("cartesian", -L/2.+np.arange(nx+1)*h, -l/2+np.arange(ny+1)*h)
else:
  m=gf.Mesh("regular_simplices", -L/2+np.arange(nx+1)*h, -l/2+np.arange(ny+1)*h)

LEFT_RG   = 101
RIGHT_RG  = 102
BOTTOM_RG = 103
TOP_RG    = 104
fleft   = m.outer_faces_with_direction([-1., 0.], 0.5)
fright  = m.outer_faces_with_direction([1., 0.], 0.5)
fbottom = m.outer_faces_with_direction([0., -1.], 0.5)
ftop    = m.outer_faces_with_direction([0., 1.], 0.5)
m.set_region(LEFT_RG,   fleft)
m.set_region(RIGHT_RG,  fright)
m.set_region(BOTTOM_RG, fbottom)
m.set_region(TOP_RG,    ftop)

m.export_to_vtu("mesh.vtu")

# Levelset definition:
R1 = 2.5
R2 = 16
ytip = R1
xtip = np.sqrt(R2*R2-R1*R1)
ls1 = gf.LevelSet(m, 2, f"y-{R1:g}*tanh(x/7.)", f"x*x+y*y-{R2*R2:g}")
ls2 = gf.LevelSet(m, 2, f"y+{R1:g}*tanh(x/7.)", f"x*x+y*y-{R2*R2:g}")
mls = gf.MeshLevelSet(m)
mls.add(ls1)
mls.add(ls2)
mls.adapt()

# Basic mesh_fem without enrichment:
mf_pre = gf.MeshFem(m)
if (quad):
  mf_pre.set_fem(gf.Fem("FEM_QK(2,2)"))
else:
  mf_pre.set_fem(gf.Fem("FEM_PK(2,2)"))

# Definition of the enriched finite element method (MeshFemLevelSet):
mfls = gf.MeshFem("levelset", mls, mf_pre)

mf_u = gf.MeshFem("clone", mfls)
#mf_u = gf.MeshFem("levelset", mls, mf_pre)
mf_u.set_qdim(2)

mf_theta = gf.MeshFem("clone", mfls)
#mf_theta = gf.MeshFem("levelset", mls, mf_pre)

# MeshIm definition (MeshImLevelSet):
if (quad):
  mim = gf.MeshIm("levelset", mls, "all",
         gf.Integ("IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(6),5)"),
         gf.Integ("IM_STRUCTURED_COMPOSITE(IM_GAUSS_PARALLELEPIPED(2,6),9)"),
         gf.Integ("IM_GAUSS_PARALLELEPIPED(2,4)"))
else:
  mim = gf.MeshIm("levelset", mls, "all",
         gf.Integ("IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(6),3)"),
         gf.Integ("IM_STRUCTURED_COMPOSITE(IM_GAUSS_PARALLELEPIPED(2,6),9)"),
         gf.Integ("IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(6),5)"))

mim_bound = [gf.MeshIm("levelset", mls, boundary,
                       gf.Integ("IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(6),3)"))
             for boundary in ["boundary(a)", "boundary(b)"]]

surf_crack = gf.asm("generic", mim_bound[0], 0, "1", -1)\
            +gf.asm("generic", mim_bound[1], 0, "1", -1)
print(f"surf_crack = {surf_crack:g}")


# Model definition:
md = gf.Model("real")
md.add_fem_variable("u", mf_u)
md.add_fem_variable("theta", mf_theta)
md.add_variable("P", 1)
md.add_variable("T", 1); md.set_variable("T", 273+room_temp)
# Data
md.add_initialized_data("surf_crack", surf_crack)
md.add_initialized_data("mu0", E0/(2*(1+nu)))
md.add_initialized_data("la0", E0*nu/((1-2*nu)*(1+nu)))
md.add_initialized_data("ka0", E0/(3*(1-2*nu)))
md.add_initialized_data("alpha", alpha)
md.add_initialized_data("alpha_env", alpha_env)
md.add_initialized_data("alpha_crack", alpha_crack)
md.add_initialized_data("beta", beta)
md.add_initialized_data("gamma", gamma)
md.add_initialized_data("mR", mR)
md.add_initialized_data("theta0", theta0)
# Lamé coefficients + bulk modulus
for p1,p2 in [["la","la0"],["mu","mu0"],["ka","ka0"]]:
  md.add_macro(p1, f"{p2}/(1+gamma*(theta-theta0))")
# Deformation gradient
md.add_macro("F", "Id(2)+Grad_u")
md.add_macro("Fm", "Id(2)+Xfem_minus(Grad_u)")
md.add_macro("Fp", "Id(2)+Xfem_plus(Grad_u)")

md.add_macro("Nanson(FF)", "Det(FF)*Norm(Inv(FF')*Normal)")
md.add_macro("th", "exp(beta*(theta-theta0))")

# Deformation tensor
md.add_macro("E", "(Grad_u+Grad_u'+Grad_u'*Grad_u)*0.5")
# Second Piola-Kirchhoff Stress tensor, elastic + thermal expansion part
md.add_macro("S", "((la*Trace(E)*Id(2)+2*mu*E)/th-1.5*(th-1/th)*ka*Id(2))")
# Elasticity term
md.add_nonlinear_term(mim, "(F*S):Grad_Test_u")
# Diffusion term
md.add_nonlinear_term(mim,
  "alpha*Det(F)*(Inv(F')*Grad_theta).(Inv(F')*Grad_Test_theta)")
for mimb in mim_bound:
  md.add_nonlinear_term(mimb, # Ideal Gas law
  "(P*((Xfem_plus(u)-Xfem_minus(u)).Normalized(Inv(Fm')*Normal)+1e-6)"
  "-mR*T/surf_crack)*Test_P")
  md.add_nonlinear_term(mimb, # Heat flux equilibrium
                        "(Nanson(Fm)*(T-273-Xfem_minus(theta))"
                        "+Nanson(Fp)*(T-273-Xfem_plus(theta)))*Test_T")
  md.add_nonlinear_term(mimb, # Heat exchange on the cracks
  "-alpha_crack*(Nanson(Fm)*(T-273-Xfem_minus(theta))*Xfem_minus(Test_theta)"
               "+Nanson(Fp)*(T-273-Xfem_plus(theta))*Xfem_plus(Test_theta))")
for var,rg in [["theta_B", BOTTOM_RG],["theta_T",TOP_RG]]: # Heat exchange at
  md.add_initialized_data(var, room_temp)                  # bottom/top
  md.add_nonlinear_term(mim,"-alpha_env*Nanson(F)*("+var+"-theta)*Test_theta",
                        rg)
for mimb in mim_bound: # Follower pressure
  md.add_nonlinear_term(mimb,"(P*Det(Fm)*Inv(Fm')*Normal).Xfem_minus(Test_u)")
  md.add_nonlinear_term(mimb,"(-P*Det(Fp)*Inv(Fp')*Normal).Xfem_plus(Test_u)")

# Fixed displacement
md.add_Dirichlet_condition_with_multipliers(mim, "u", 1, LEFT_RG)
md.add_Dirichlet_condition_with_multipliers(mim, "u", 1, RIGHT_RG)


# Increase temperature gap
for it,fact in enumerate([0,1.,2.,3.,3.5,4.]):
#  md.set_variable("theta_T", room_temp + fact*jump_temp)
  md.set_variable("theta_T", room_temp + 0.2*jump_temp)
  if it > 0:
    ls_new = gf.LevelSet(m, 2, f"y-{R1:g}*tanh(x/7.)", f"x*x+y*y-{(R2+fact)**2:g}")
    ls1.set_values(ls_new.values(0),ls_new.values(1))
    mls.adapt()
    mf_u.adapt()
    mf_theta.adapt()
    mim.adapt()
    md.set_variable("u", 0*md.variable("u"))
    md.set_variable("theta", 0*md.variable("theta"))
  print(f"Displacement dofs: {mf_u.nbdof()}\nTotal model dofs: {md.nbdof()}")

  # Solve with fixed pressure (to open the crack lips)
  md.disable_variable("P");
  md.set_variable("P", [1e-2])
  md.solve("max_res", 5E-4, "max_iter", 100, "noisy")
  # Solve fully coupled problem
  md.enable_variable("P");

  md.solve("max_res", 5E-8, "max_iter", 100, "noisy")

  U = md.variable("u")
  theta = md.variable("theta")
  T = md.variable("T")[0]
  print(f"Gas temperature {T-273:g} C")
  P = md.variable("P")[0]
  print(f"Gas pressure {P:g} MPa")

  # Interpolation of the solution on a cut mesh for the drawing purpose
  cut_mesh = mls.cut_mesh()
  mfv = gf.MeshFem(cut_mesh, 2)
  mfv.set_classical_discontinuous_fem(2, 0.001)
  mfw = gf.MeshFem(cut_mesh, 1)
  mfw.set_classical_discontinuous_fem(2, 0.001)

  # For the computation of the Von Mises stress
  mfvm = gf.MeshFem(cut_mesh)
  mfvm.set_classical_discontinuous_fem(4, 0.001)
  md.add_interpolate_transformation_from_expression(f"id_trans{it}", cut_mesh, m, "X")
  md.add_macro(f"gradu_{it}", f"Interpolate(Grad_u, id_trans{it})")
  md.add_macro(f"theta_{it}", f"Interpolate(theta, id_trans{it})")
  md.add_macro(f"F_{it}", f"Id(2)+gradu_{it}")
  md.add_macro(f"E_{it}", f"gradu_{it}+gradu_{it}'+gradu_{it}'*gradu_{it}")
  md.add_macro(f"mu_{it}", f"mu0/(1+gamma*(theta_{it}-theta0))")
  md.add_macro(f"la_{it}", f"la0/(1+gamma*(theta_{it}-theta0))")
  md.add_macro(f"ka_{it}", f"ka0/(1+gamma*(theta_{it}-theta0))")
  md.add_macro(f"th_{it}", f"exp(beta*(theta_{it}-theta0))")
  md.add_macro(f"S_{it}", f"((la_{it}*Trace(E_{it})*Id(2)+2*mu_{it}*E_{it})/th_{it}-1.5*(th_{it}-1/th_{it})*ka_{it}*Id(2))")

  V  = gf.compute_interpolate_on(mf_u, U, mfv)
  Th = gf.compute_interpolate_on(mf_theta, theta, mfw)
  VM = md.interpolation(f"sqrt(3/2)*Norm(Deviator(Cauchy_stress_from_PK2(S_{it}, gradu_{it})))", mfvm)

  mfv.export_to_vtu(f"cracked_body_{it}.vtu", V, "V", mfw, Th, "Temperature", mfvm, VM, "Von Mises stress")

print("You can view the solution with (for example):")
print("paraview cracked_body_2.vtu")
1 Like
4

Notes mainly to myself …

What was a bit unclear to me for supporting the use case described in Janine’s original post is the adaptation of the mf_sing_u mesh_fem. This object is created with

 mf_sing_u = gf.MeshFem(‘global function’,m,ls,[ckoff0,ckoff1,ckoff2,ckoff3],1)

It contains 4 global functions defined in the coordinate system described by the levelset ls. When ls is adapted, these functions do not need to be adapted as well. Every time these functions are called they use the current version of the ls levelset object they are linked to.

Internally mf_sing_u is a mesh_fem_global_function C++ object, which allocates and owns a fem_global_function object, which in turn contains a link to the 4 base functions in a vector of pointers std::shared_ptr<const global_function>.

The mesh_fem_global_function object is a thin object that doesn’t seem to need updating.

The fem_global_function object stored in the mesh_fem_global_function object, is a thick object containing all the information about which base functions are active in which elements. Hence, when the base functions change, because the underlying levelset has changed, the fem_global_function object using them, must be adapted. This object has an update_from_context function that should do the job. But I am not sure if it is triggered

The function update_from_context is called in GetFEM normally in the object constructor and then inside the go_check function of all context aware objects, which in turn is called from the context_check function, if an object context state has been marked as other than normal.

For the fem_global_function object specifically context_check is called whenever nb_dof or index_of_global_dof is queried. So, provided that the context state of the object has been marked as not normal, the update of this object should occur automagically next time this object is used.

The only problem now is that the fem_global_function object is not aware of its dependence on the levelset object that its base functions depend on. This is because it links to the base class global_function which is not context aware and not the specific implementation global_function_on_levelsets_2D_ which is in fact context aware. Hence the information that levelset ls has changed will not propagate all the way to the fem_global_function object.

1 Like
5

This issue must now be fixed with the latest commits [1,2] in the master branch.

@jab2443 could you test it?

Otherwise, I have tested this version with the following script

import sys
import numpy as np
import getfem as gf

np.set_printoptions(threshold=sys.maxsize)
gf.util_trace_level(1)

# Parameters:  #########################################################
E0 = 5E2            # [N/mm^2] Young modulus at reference temperature
nu = 0.45           # [-] Poisson ratio
beta = 8e-5         # [1/K] Expansion coefficient
gamma = 0.1         # [1/K] Coefficient of thermal weakening
theta0 = 30.        # [C] Reference temperature
alpha = 0.1         # [mm^2/s] Thermal diffusivity (0.1 for rubber)
alpha_env = 0.02    # [N/mm/s/K] Heat transfer coefficient
alpha_crack = 0.01  # [N/mm/s/K] Heat transfer coefficient
mR = 287e3*5e-9     # [N*mm/K] --> 287000 [N*mm/kg/K] * mass [kg]
room_temp = 30.     # [C] temperature on one side of the block
jump_temp = 100.    # [C] temperature jump between the two sides

# Mesh definition:
L = 40.             # [mm]
l = 10.             # [mm]
ny = 10             # Number of elements in each direction
quad = True         # quad mesh or triangle one

h = l/ny
nx = np.floor(L/h)

if (quad):
  m=gf.Mesh("cartesian", -L/2.+np.arange(nx+1)*h, -l/2+np.arange(ny+1)*h)
else:
  m=gf.Mesh("regular_simplices", -L/2+np.arange(nx+1)*h, -l/2+np.arange(ny+1)*h)

LEFT_RG   = 101
RIGHT_RG  = 102
BOTTOM_RG = 103
TOP_RG    = 104
fleft   = m.outer_faces_with_direction([-1., 0.], 0.5)
fright  = m.outer_faces_with_direction([1., 0.], 0.5)
fbottom = m.outer_faces_with_direction([0., -1.], 0.5)
ftop    = m.outer_faces_with_direction([0., 1.], 0.5)
m.set_region(LEFT_RG,   fleft)
m.set_region(RIGHT_RG,  fright)
m.set_region(BOTTOM_RG, fbottom)
m.set_region(TOP_RG,    ftop)

m.export_to_vtu("mesh.vtu")

# Levelset definition:
R1 = 2.5
R2 = 16.
# approximate tip position
ytip = R1
xtip = np.sqrt(R2*R2-R1*R1)
ls1 = gf.LevelSet(m, 2, f"y-{R1:g}*tanh(x/7.)", f"x*x+y*y-{R2*R2:g}")
ls2 = gf.LevelSet(m, 2, f"y+{R1:g}*tanh(x/7.)", f"x*x+y*y-{R2*R2:g}")
mls = gf.MeshLevelSet(m)
mls.add(ls1)
mls.add(ls2)
mls.adapt()

# Basic mesh_fem without enrichment:
mf_pre = gf.MeshFem(m)
if (quad):
  mf_pre.set_fem(gf.Fem("FEM_QK(2,2)"))
else:
  mf_pre.set_fem(gf.Fem("FEM_PK(2,2)"))

# Definition of the enriched finite element method (MeshFemLevelSet):
mfls = gf.MeshFem("levelset", mls, mf_pre)

# Global functions for asymptotic enrichment:
mf_PoU = gf.MeshFem(m)
mf_PoU.set_classical_fem(1) # partition of unity
ck = [gf.GlobalFunction("crack",i) for i in range(4)]
mf_enr = [gf.MeshFem("global function", m, ls, ck, 1) for ls in [ls1,ls2]]
mf_sing = [gf.MeshFem("product", mf_PoU, mf) for mf in mf_enr]
if True:
  DOFpts = mf_PoU.basic_dof_nodes()
  ctip_dofs = [np.nonzero(np.linalg.norm(DOFpts-x,axis=0) < 1.)[0]
               for x in [[[xtip],[-ytip]],[[-xtip],[ytip] ],
                         [[xtip],[ytip]], [[-xtip],[-ytip]]]]
  mf_sing[0].set_enriched_dofs(np.union1d(ctip_dofs[0],ctip_dofs[1]))
  mf_sing[1].set_enriched_dofs(np.union1d(ctip_dofs[2],ctip_dofs[3]))
mf_u = gf.MeshFem("sum", mf_sing[0], mf_sing[1], mfls)
mf_u.set_qdim(2)

mf_theta = gf.MeshFem("sum", mf_sing[0], mf_sing[1], mfls)

# MeshIm definition (MeshImLevelSet):
if (quad):
  mim = gf.MeshIm("levelset", mls, "all",
         gf.Integ("IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(6),5)"),
         gf.Integ("IM_STRUCTURED_COMPOSITE(IM_GAUSS_PARALLELEPIPED(2,6),9)"),
         gf.Integ("IM_GAUSS_PARALLELEPIPED(2,4)"))
else:
  mim = gf.MeshIm("levelset", mls, "all",
         gf.Integ("IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(6),3)"),
         gf.Integ("IM_STRUCTURED_COMPOSITE(IM_GAUSS_PARALLELEPIPED(2,6),9)"),
         gf.Integ("IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(6),5)"))

mim_bound = [gf.MeshIm("levelset", mls, boundary,
                       gf.Integ("IM_STRUCTURED_COMPOSITE(IM_TRIANGLE(6),3)"))
             for boundary in ["boundary(a)", "boundary(b)"]]

surf_crack = gf.asm("generic", mim_bound[0], 0, "1", -1)\
            +gf.asm("generic", mim_bound[1], 0, "1", -1)
print(f"surf_crack = {surf_crack:g}")


# Model definition:
md = gf.Model("real")
md.add_fem_variable("u", mf_u)
md.add_fem_variable("theta", mf_theta)
md.add_variable("P", 1)
md.add_variable("T", 1)
md.set_variable("T", 273+room_temp)
# Data
md.add_initialized_data("surf_crack", surf_crack)
md.add_initialized_data("mu0", E0/(2*(1+nu)))
md.add_initialized_data("la0", E0*nu/((1-2*nu)*(1+nu)))
md.add_initialized_data("ka0", E0/(3*(1-2*nu)))
md.add_initialized_data("alpha", alpha)
md.add_initialized_data("alpha_env", alpha_env)
md.add_initialized_data("alpha_crack", alpha_crack)
md.add_initialized_data("beta", beta)
md.add_initialized_data("gamma", gamma)
md.add_initialized_data("mR", mR)
md.add_initialized_data("theta0", theta0)
# Lamé coefficients + bulk modulus
for p1,p2 in [["la","la0"],["mu","mu0"],["ka","ka0"]]:
  md.add_macro(p1, f"{p2}/(1+gamma*(theta-theta0))")
# Deformation gradient
md.add_macro("F", "Id(2)+Grad_u")
md.add_macro("Fm", "Id(2)+Xfem_minus(Grad_u)")
md.add_macro("Fp", "Id(2)+Xfem_plus(Grad_u)")

md.add_macro("Nanson(FF)", "Det(FF)*Norm(Inv(FF')*Normal)")
md.add_macro("th", "exp(beta*(theta-theta0))")

# Deformation tensor
md.add_macro("E", "(Grad_u+Grad_u'+Grad_u'*Grad_u)*0.5")
# Second Piola-Kirchhoff Stress tensor, elastic + thermal expansion part
md.add_macro("S", "((la*Trace(E)*Id(2)+2*mu*E)/th-1.5*(th-1/th)*ka*Id(2))")
# Elasticity term
md.add_nonlinear_term(mim, "(F*S):Grad_Test_u")
# Diffusion term
md.add_nonlinear_term(mim,
  "alpha*Det(F)*(Inv(F')*Grad_theta).(Inv(F')*Grad_Test_theta)")
for mimb in mim_bound:
  md.add_nonlinear_term(mimb, # Ideal Gas law
  "(P*((Xfem_plus(u)-Xfem_minus(u)).Normalized(Inv(Fm')*Normal)+1e-6)"
  "-mR*T/surf_crack)*Test_P")
  md.add_nonlinear_term(mimb, # Heat flux equilibrium
                        "(Nanson(Fm)*(T-273-Xfem_minus(theta))"
                        "+Nanson(Fp)*(T-273-Xfem_plus(theta)))*Test_T")
  md.add_nonlinear_term(mimb, # Heat exchange on the cracks
  "-alpha_crack*(Nanson(Fm)*(T-273-Xfem_minus(theta))*Xfem_minus(Test_theta)"
               "+Nanson(Fp)*(T-273-Xfem_plus(theta))*Xfem_plus(Test_theta))")
for var,rg in [["theta_B", BOTTOM_RG],["theta_T",TOP_RG]]: # Heat exchange at
  md.add_initialized_data(var, room_temp)                  # bottom/top
  md.add_nonlinear_term(mim,"-alpha_env*Nanson(F)*("+var+"-theta)*Test_theta",
                        rg)
for mimb in mim_bound: # Follower pressure
  md.add_nonlinear_term(mimb,"(P*Det(Fm)*Inv(Fm')*Normal).Xfem_minus(Test_u)")
  md.add_nonlinear_term(mimb,"(-P*Det(Fp)*Inv(Fp')*Normal).Xfem_plus(Test_u)")

# Fixed displacement
md.add_Dirichlet_condition_with_multipliers(mim, "u", 1, LEFT_RG)
md.add_Dirichlet_condition_with_multipliers(mim, "u", 1, RIGHT_RG)


md.set_variable("theta_T", room_temp + 0.2*jump_temp)
for it,fact in enumerate([0,1.,2.,3.,3.25,3.5]):
  if it > 0:
    #mf_u = mf_sing[0]       + mf_sing[1]       + mfls
    #       mf_PoU*mf_enr[0] + mf_PoU*mf_enr[1] + mfls(mls,mf_pre)

    # approximate tip position
    xtip = np.sqrt((R2+fact)**2-R1**2)
    ls_new = gf.LevelSet(m, 2, f"y-{R1:g}*tanh(x/7.)", f"x*x+y*y-{(R2+fact)**2:g}")
    ls1.set_values(ls_new.values(0),ls_new.values(1))
    mls.adapt()
    mfls.adapt()
    ctip_dofs = [np.nonzero(np.linalg.norm(DOFpts-x,axis=0) < 1.)[0]
                 for x in [[[xtip],[-ytip]],[[-xtip],[ytip] ]]]
    mf_sing[0].set_enriched_dofs(np.union1d(ctip_dofs[0],ctip_dofs[1]))
    mf_sing[0].adapt()
    mf_u.adapt()
    mf_theta.adapt()
    mim.adapt()
    md.set_variable("u", 0*md.variable("u"))
    md.set_variable("theta", 0*md.variable("theta"))
  print(f"Displacement dofs: {mf_u.nbdof()}\nTotal model dofs: {md.nbdof()}")

  # Solve with fixed pressure (to open the crack lips)
  md.disable_variable("P")
  md.set_variable("P", [1e-2])
  md.solve("max_res", 5E-4, "max_iter", 100, "noisy")
  # Solve fully coupled problem
  md.enable_variable("P")

  md.solve("max_res", 5E-8, "max_iter", 100, "noisy")

  U = md.variable("u")
  theta = md.variable("theta")
  T = md.variable("T")[0]
  print(f"Gas temperature {T-273:g} C")
  P = md.variable("P")[0]
  print(f"Gas pressure {P:g} MPa")

  # Interpolation of the solution on a cut mesh for the drawing purpose
  cut_mesh = mls.cut_mesh()
  mfv = gf.MeshFem(cut_mesh, 2)
  mfv.set_classical_discontinuous_fem(2, 0.001)
  mfw = gf.MeshFem(cut_mesh, 1)
  mfw.set_classical_discontinuous_fem(2, 0.001)

  # For the computation of the Von Mises stress
  mfvm = gf.MeshFem(cut_mesh)
  mfvm.set_classical_discontinuous_fem(4, 0.001)
  md.add_interpolate_transformation_from_expression(f"id_trans{it}", cut_mesh, m, "X")
  md.add_macro(f"gradu_{it}", f"Interpolate(Grad_u, id_trans{it})")
  md.add_macro(f"theta_{it}", f"Interpolate(theta, id_trans{it})")
  md.add_macro(f"F_{it}", f"Id(2)+gradu_{it}")
  md.add_macro(f"E_{it}", f"gradu_{it}+gradu_{it}'+gradu_{it}'*gradu_{it}")
  md.add_macro(f"mu_{it}", f"mu0/(1+gamma*(theta_{it}-theta0))")
  md.add_macro(f"la_{it}", f"la0/(1+gamma*(theta_{it}-theta0))")
  md.add_macro(f"ka_{it}", f"ka0/(1+gamma*(theta_{it}-theta0))")
  md.add_macro(f"th_{it}", f"exp(beta*(theta_{it}-theta0))")
  md.add_macro(f"S_{it}", f"((la_{it}*Trace(E_{it})*Id(2)+2*mu_{it}*E_{it})/th_{it}-1.5*(th_{it}-1/th_{it})*ka_{it}*Id(2))")

  V  = gf.compute_interpolate_on(mf_u, U, mfv)
  Th = gf.compute_interpolate_on(mf_theta, theta, mfw)
  VM = md.interpolation(f"sqrt(3/2)*Norm(Deviator(Cauchy_stress_from_PK2(S_{it}, gradu_{it})))", mfvm)

  mfv.export_to_vtu(f"cracked_body_{it}.vtu", V, "V", mfw, Th, "Temperature", mfvm, VM, "Von Mises stress")

print("You can view the solution with (for example):")
print("paraview cracked_body_2.vtu")
1 Like