Adapt MeshFem with XFEM enrichment

[Janine Ashleigh Birnbaum]

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 = ck0*coff # gf.GlobalFunction('product',ck0,coff)
ckoff1 = ck1*coff
ckoff2 = ck2*coff
ckoff3 = ck3*coff

# 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+2*Mu/(Lambda+2*Mu); 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 /= 2*np.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()