```
import os
import numpy as np
import matplotlib.pyplot as plt
import dolfin as dol
import meshio
def convert_msh_to_XDMF(filename, dim, dim_keep = None):
'''reads a .msh file and writes mesh, subdomains, and boundary subdomains
to XDMF files'''
if filename.split(sep = '.')[-1] != 'msh':
raise TypeError('.msh file required')
fn = '.'.join(filename.split(sep ='.')[:-1])
if dim_keep is None:
dim_keep = dim
msh = meshio.read(filename)
if dim == 1:
volume_cell_type = 'line'
surface_cell_type = 'vertex'
elif dim == 2:
volume_cell_type = 'triangle'
surface_cell_type = 'line'
elif dim == 3:
volume_cell_type = 'tetra'
surface_cell_type = 'triangle'
else:
raise ValueError('dim can only be 1,2, or 3')
volume_cells = []
surface_cells = []
for cell in msh.cells:
if cell.type == volume_cell_type:
if len(volume_cells) == 0:
volume_cells = cell.data
else:
volume_cells = np.vstack([volume_cells, cell.data])
if cell.type == surface_cell_type:
if len(surface_cells) == 0:
surface_cells = cell.data
else:
surface_cells = np.vstack([surface_cells, cell.data])
volume_data = None
surface_data = None
for key in msh.cell_data_dict['gmsh:physical'].keys():
if key == volume_cell_type:
volume_data = msh.cell_data_dict['gmsh:physical'][key]
elif key == surface_cell_type:
surface_data = msh.cell_data_dict['gmsh:physical'][key]
volume_mesh = meshio.Mesh(points = msh.points[:,:dim_keep],
cells = [(volume_cell_type, volume_cells)],
cell_data = {'marker' : [volume_data]})
meshio.write(fn + '_{}D_mesh.xdmf'.format(dim), volume_mesh)
if not (surface_data is None):
surface_mesh = meshio.Mesh(points = msh.points[:,:dim_keep],
cells = [(surface_cell_type, surface_cells)],
cell_data = {'marker' : [surface_data]})
meshio.write(fn + '_{}D_mesh.xdmf'.format(dim-1), surface_mesh)
def load_XDMF_files(filename, dim):
'''create mesh from XDMF files'''
mesh = dol.Mesh()
fn = '.'.join(filename.split(sep ='.')[:-1])
with dol.XDMFFile(fn + '_{}D_mesh.xdmf'.format(dim)) as infile:
infile.read(mesh)
mvc = dol.MeshValueCollection('size_t', mesh = mesh, dim = dim)
with dol.XDMFFile(fn + '_{}D_mesh.xdmf'.format(dim)) as infile:
infile.read(mvc, 'marker')
volume_marker = dol.cpp.mesh.MeshFunctionSizet(mesh, mvc)
del(mvc)
dx = dol.Measure('dx', domain = mesh, subdomain_data = volume_marker)
if os.path.exists(fn + '_{}D_mesh.xdmf'.format(dim-1)):
mvc = dol.MeshValueCollection('size_t', mesh = mesh, dim = dim -1)
with dol.XDMFFile(fn + '_{}D_mesh.xdmf'.format(dim-1)) as infile:
infile.read(mvc, 'marker')
boundary_marker = dol.cpp.mesh.MeshFunctionSizet(mesh, mvc)
del(mvc)
ds = dol.Measure('ds', domain = mesh, subdomain_data = boundary_marker)
dS = dol.Measure('dS', domain = mesh, subdomain_data = boundary_marker)
else:
warn('{}D mesh not found'.format(dim-1))
boundary_marker = None
ds = dol.Measure('ds', domain = mesh)
dS = dol.Measure('dS', domain = mesh)
return {'mesh' : mesh,
'volume_measure' : dx,
'boundary_measure' : ds,
'internal_boundary_measure' : dS,
'volume_marker' : volume_marker,
'boundary_marker' : boundary_marker}
class Rho(dol.UserExpression):
def __init__(self, phi):
self.phi = phi
dol.UserExpression.__init__(self)
def eval(self, value, x):
print(x[2]) # I added this to check where the function is evaluated
mu = self.phi(x)
value[0] = - mu**2 * np.sign(mu)
def value_shape(self):
return ()
class Rho_Der(dol.UserExpression):
def __init__(self, phi):
self.phi = phi
dol.UserExpression.__init__(self)
def eval(self, value, x):
mu = self.phi(x)
value[0] = - 2. * abs(mu)
def value_shape(self):
return ()
#Actual program starts here
convert_msh_to_XDMF( 'structure.msh', dim = 3)
mesh_dic = load_XDMF_files('structure.msh', dim = 3)
V = dol.FunctionSpace(
mesh_dic['mesh'],
"P", 1, #linear Lagrange elements
)
# trial and test functions
u = dol.Function(V)
v = dol.TestFunction(V)
du = dol.TrialFunction(V)
# measures
dx = mesh_dic['volume_measure']
ds = mesh_dic['boundary_measure']
dS = mesh_dic['internal_boundary_measure']
# quadratic form and its jacobian
a = dol.inner(dol.grad(u), dol.grad(v)) * dx
da = dol.inner(dol.grad(du), dol.grad(v)) * dx
nl = dol.avg(v) * Rho(u) * dS(3)
dnl = dol.avg(v) * dol.avg(du) * Rho_Der(u) * dS(3)
bc = []
bc.append(
dol.DirichletBC(V,
dol.Constant(0),
mesh_dic['boundary_marker'],
1))
bc.append(
dol.DirichletBC(V,
dol.Constant(1),
mesh_dic['boundary_marker'],
2)
)
# solve non-linear problem
functional = a - nl
jacobian = da - dnl
problem = dol.NonlinearVariationalProblem(functional, u, bc, jacobian)
solver = dol.NonlinearVariationalSolver(problem)
solver.solve()
```

Here is the code but to run it you will need the .msh file.

How can I attach it?

Best.

Iacopo