Hello everyone, I’ve been working with the MeshTie demos for thermo-mechanical problems (asimov-contact/python/demos/meshtie_demos/thermo_mech.py at main · Wells-Group/asimov-contact · GitHub), and running into an issue when I substitute a custom gmsh mesh. The setup is identical to the demo except for the mesh, I am using a sheet with an electrode above, coupled across an interface using MeshTie. With an applied temperature of 1, I would expect the solution to stay within [0, 1], but I am seeing lower values going negative and higher up to ~1.06 from the very first timestep.
The code is (For simplification I have considered only thermal part)
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
from dolfinx import default_scalar_type, log
from dolfinx.fem import ( Function,form, functionspace)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix,
assemble_vector, create_vector, set_bc)
from ufl import (Measure, TestFunction, TrialFunction,
grad, inner, lhs, rhs)
from dolfinx.graph import adjacencylist
from dolfinx.io import VTXWriter
from dolfinx_contact.cpp import MeshTie, Problem
from dolfinx_contact.parallel_mesh_ghosting import create_contact_mesh
import gmsh
import dolfinx.io.gmsh
import dolfinx.fem as _fem
def create_electrode_sheet_mesh(
model,
res=0.01,
comm: MPI.Comm = MPI.COMM_WORLD,
rank: int = 0,
):
if comm.rank == rank:
gmsh.option.setNumber("General.Terminal", 0)
gmsh.option.setNumber("Mesh.CharacteristicLengthFactor", res)
# Sheet rectangle: (0,0)–(50,3)
p1 = model.occ.addPoint(0.0, 0.0, 0.0)
p2 = model.occ.addPoint(50.0, 0.0, 0.0)
p3 = model.occ.addPoint(50.0, 3.0, 0.0)
p4 = model.occ.addPoint(28.0, 3.0, 0.0)
p5 = model.occ.addPoint(22.0, 3.0, 0.0)
p6 = model.occ.addPoint(0.0, 3.0, 0.0)
# Electrode polygon
p7 = model.occ.addPoint(20.0, 16.0, 0.0)
p8 = model.occ.addPoint(30.0, 16.0, 0.0)
p9 = model.occ.addPoint(20.0, 6.0, 0.0)
p10 = model.occ.addPoint(30.0, 6.0, 0.0)
p11 = model.occ.addPoint(22.0, 3.01, 0.0)
p12 = model.occ.addPoint(28.0, 3.01, 0.0)
l6 = model.occ.addLine(p5, p6) # 5 -> 6
l1 = model.occ.addLine(p6, p1) # 6 -> 1
l2 = model.occ.addLine(p1, p2) # 1 -> 2
l3 = model.occ.addLine(p2, p3) # 2 -> 3
l4 = model.occ.addLine(p3, p4) # 3 -> 4
l5 = model.occ.addLine(p4, p5) # 4 -> 5
l7 = model.occ.addLine(p12, p11) # 12 -> 11
l8 = model.occ.addLine(p11, p9) # 11 -> 9
l9 = model.occ.addLine(p9, p7) # 9 -> 7
l10 = model.occ.addLine(p7, p8) # 7 -> 8
l11 = model.occ.addLine(p8, p10) # 8 -> 10
l12 = model.occ.addLine(p10, p12) # 10 -> 12
cl1 = model.occ.addCurveLoop([l7, l8, l9, l10, l11, l12])
s1 = model.occ.addPlaneSurface([cl1]) # Electrode
cl2 = model.occ.addCurveLoop([l6, l1, l2, l3, l4, l5])
s2 = model.occ.addPlaneSurface([cl2]) # Sheet
model.occ.synchronize()
# Physical groups (like your .geo)
model.addPhysicalGroup(2, [s1], tag=13)
model.setPhysicalName(2, 13, "Electrode")
model.addPhysicalGroup(2, [s2], tag=14)
model.setPhysicalName(2, 14, "Sheet")
# Curves
model.addPhysicalGroup(1, [l7], tag=15)
model.setPhysicalName(1, 15, "Contact_surface_electrode")
model.addPhysicalGroup(1, [l5], tag=16)
model.setPhysicalName(1, 16, "Contact_surface_sheet")
model.addPhysicalGroup(1, [l2], tag=17)
model.setPhysicalName(1, 17, "Fixed_bottom")
model.addPhysicalGroup(1, [l10], tag=18)
model.setPhysicalName(1, 18, "Loading_surface")
model.mesh.generate(2)
return model
else:
return None
#Creating mesh from gmsh
gmsh.initialize()
name = "sheet_and_electrode"
model = gmsh.model()
model.add(name)
model.setCurrent(name)
model = create_electrode_sheet_mesh(model, res=0.01)
mesh_data = dolfinx.io.gmsh.model_to_mesh(model, MPI.COMM_WORLD, 0, gdim=2)
gmsh.finalize()
mesh = mesh_data.mesh
# read in mesh and markers
tdim = mesh.topology.dim
gdim = mesh.geometry.dim
mesh.topology.create_connectivity(tdim - 1, 0)
mesh.topology.create_connectivity(tdim - 1, tdim)
domain_marker = mesh_data.cell_tags
facet_marker = mesh_data.facet_tags
# contact_bdy_1 = 15 #Contact surface of electrode
# contact_bdy_2 = 16 #Contact surface of Sheet
if mesh.comm.size > 1:
mesh, facet_marker, domain_marker = create_contact_mesh(
mesh, facet_marker, domain_marker, [15,16]
)
# compiler options to improve performance
jit_options = {"cffi_extra_compile_args": [], "cffi_libraries": ["m"]}
# Linear solver options
ksp_tol = 1e-10
petsc_options = {
"ksp_type": "cg",
"ksp_rtol": ksp_tol,
"ksp_atol": ksp_tol,
"pc_type": "gamg",
"mg_levels_ksp_type": "chebyshev",
"mg_levels_pc_type": "jacobi",
"pc_gamg_type": "agg",
"pc_gamg_coarse_eq_limit": 100,
"pc_gamg_agg_nsmooths": 1,
"pc_gamg_threshold": 1e-3,
"pc_gamg_square_graph": 2,
"ksp_norm_type": "unpreconditioned",
}
# measures
dx = Measure("dx", domain=mesh, subdomain_data=domain_marker)
ds = Measure("ds", domain=mesh, subdomain_data=facet_marker)
# Thermal problem ufl (implicit Euler)
Q = functionspace(mesh, ("CG", 1))
q, r = TrialFunction(Q), TestFunction(Q)
T0 = Function(Q)
kdt = 0.01
Q0 = functionspace(mesh, ("DG", 0))
#MeshTie requires ALL coefficients to be Functions, not scalars so kdt_custom as function
kdt_custom = Function(Q0)
kdt_custom.interpolate(lambda x: np.full((1, x.shape[1]), kdt))
#Weak form of thermal equation
therm = (q - T0) * r * dx + kdt * inner(grad(q), grad(r)) * dx
a_therm, L_therm = lhs(therm), rhs(therm)
T0.x.array[:] = 0
#Setting boundary conditions for thermal problem at bottom surface (T=1 at y=0)
dofs = _fem.locate_dofs_topological(Q, entity_dim=tdim - 1, entities=facet_marker.find(17)) # Bottom line
Tbc = _fem.dirichletbc(value=default_scalar_type((0.0)), dofs=dofs, V=Q)
dofs2 = _fem.locate_dofs_topological(Q, entity_dim=tdim - 1, entities=facet_marker.find(15)) #Loading line
Tbc2 = _fem.dirichletbc(value=default_scalar_type((1.0)), dofs=dofs2, V=Q)
# surface data for Nitsche
gamma = 10
theta = 1
contact = [(0, 1), (1, 0)]
data = np.array([15, 16], dtype=np.int32)
offsets = np.array([0, 2], dtype=np.int32)
surfaces = adjacencylist(data, offsets)
# initialise meshties
meshties = MeshTie(
[facet_marker._cpp_object], surfaces._cpp_object, contact, mesh._cpp_object, quadrature_degree=3
)
meshties.generate_kernel_data(
Problem.Poisson,
Q._cpp_object,
{"T": T0._cpp_object, "kdt": kdt_custom._cpp_object},
gamma,
theta,
)
# Create matrix and vector
a_therm = form(a_therm, jit_options=jit_options)
L_therm = form(L_therm, jit_options=jit_options)
mat_therm = meshties.create_matrix(a_therm._cpp_object)
#vec_therm = create_vector(L_therm)
V = L_therm.function_spaces[0]
vec_therm = create_vector(V)
# Thermal problem: functions for updating matrix and vector
def assemble_mat_therm(A):
A.zeroEntries()
meshties.assemble_matrix(A, Q._cpp_object, Problem.Poisson)
assemble_matrix(A, a_therm, bcs=[Tbc, Tbc2])
A.assemble()
def assemble_vec_therm(b):
b.zeroEntries()
b.ghostUpdate(
addv=PETSc.InsertMode.INSERT, # type: ignore
mode=PETSc.ScatterMode.FORWARD,
) # type: ignore
assemble_vector(b, L_therm)
# Apply boundary condition and scatter reverse
apply_lifting(b, [a_therm], bcs=[[Tbc, Tbc2]], alpha=1.0)
b.ghostUpdate(
addv=PETSc.InsertMode.ADD, # type: ignore
mode=PETSc.ScatterMode.REVERSE,
) # type: ignore
set_bc(b, [Tbc, Tbc2])
# Thermal problem: create linear solver
ksp_therm = PETSc.KSP().create(mesh.comm) # type: ignore
prefix_therm = "Solver_thermal_"
ksp_therm.setOptionsPrefix(prefix_therm)
opts = PETSc.Options() # type: ignore
opts.prefixPush(ksp_therm.getOptionsPrefix())
for key in petsc_options:
opts[key] = petsc_options[key]
opts.prefixPop()
ksp_therm.setFromOptions()
ksp_therm.setOperators(mat_therm)
ksp_therm.setMonitor(lambda _, its, rnorm: print(f"Iteration: {its}, rel. residual: {rnorm}"))
time_steps = 5
T0.name = "temperature"
vtx = VTXWriter(mesh.comm, "results/thermo_meshtie.bp", [T0], "bp4")
vtx.write(0)
for i in range(time_steps):
assemble_mat_therm(mat_therm)
assemble_vec_therm(vec_therm)
ksp_therm.solve(vec_therm, T0.x.petsc_vec)
T0.x.scatter_forward()
T_array = T0.x.array
print(f"Step {i+1}: min={T_array.min():.6f}, max={T_array.max():.6f}")
vtx.write(i + 1)
vtx.close()
The demo works correctly with its built-in mesh. Has anyone encountered this with custom meshes, or is there a recommended way to look for MeshTie problems? Any suggestion is appreciated. Thanks for your time.