Dear dokken,
I tried the method you provided and plotted the normal vector on the interface. I found an issue, but I don’t know what the issue is. Is there a simple way to restrict which direction the normal vector points to, like in the old dolfin:
a = inner(u, v)('+')*dS + + Constant(0)*inner(u,v)*dx
I saw that you provided a one-sided integral entities solution to circumvent this issue, but I will need to project the jump qunatities u(“+”) - u(“-”) on the normal, so I am not sure if this solution could be used here.
Despite the arbitrary direction, I found that the magnitude is incorrect for a few facets in the adaptation of your previous solution. I’m not sure what went wrong with the ellipse. I tried a circle and got the correct result. I tried the old dolfin version, and the results were correct in both cases. Could you help me with this problem or give me a direction? Thank you! Here is my code:
import time
import dolfinx.fem.petsc
import numpy as np
import ufl
from dolfinx import cpp, fem
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.io import XDMFFile
from ufl import Measure
import gmsh
from dolfinx.io import gmshio
gmsh.initialize()
gmsh.model.add("square_with_circle")
occ = gmsh.model.occ
# --- Parameters ---
h = 20
h2 = 0.5 # 1
L_half = 500
minor_axis = 40 # 4
major_axis = 80 # 8
# Ellipse orientation (90° domain)
minor_axis_x = minor_axis
minor_axis_y = 0
major_axis_x = 0
major_axis_y = major_axis
# --- Ellipse-defining points ---
occ.addPoint(minor_axis_x, -minor_axis_y, 0, h2, tag=1)
occ.addPoint(-minor_axis_x, minor_axis_y, 0, h2, tag=2)
occ.addPoint(major_axis_x, major_axis_y, 0, h2, tag=3)
occ.addPoint(-major_axis_x, -major_axis_y, 0, h2, tag=4)
# --- Square corners ---
occ.addPoint(L_half, L_half, 0, h, tag=5)
occ.addPoint(-L_half, L_half, 0, h, tag=6)
occ.addPoint(-L_half, -L_half, 0, h, tag=7)
occ.addPoint(L_half, -L_half, 0, h, tag=8)
# --- Square edges ---
occ.addLine(6, 5, tag=1) # top
occ.addLine(5, 8, tag=2) # right
occ.addLine(8, 7, tag=3) # bottom
occ.addLine(7, 6, tag=4) # left
# --- Center of ellipse ---
occ.addPoint(0, 0, 0, h, tag=9)
# --- Elliptical arcs ---
occ.addEllipseArc(1, 9, 3, 3, tag=5) # (p1->p3 via center)
occ.addEllipseArc(3, 9, 3, 2, tag=6) # (p3->p2)
occ.addEllipseArc(2, 9, 4, 4, tag=7) # (p2->p4)
occ.addEllipseArc(4, 9, 4, 1, tag=8) # (p4->p1)
# --- Curve loops ---
outer_loop = occ.addCurveLoop([1, 2, 3, 4], tag=1) # square
ellipse_loop = occ.addCurveLoop([6, 7, 8, 5], tag=2) # ellipse (consistent order)
# --- Surfaces ---
# Matrix = square with ellipse inside (like Plane Surface(1) = {1,2,3})
s_matrix = occ.addPlaneSurface([outer_loop, ellipse_loop], tag=1)
# Inclusion = ellipse alone (like Plane Surface(2) = {4})
s_incl = occ.addPlaneSurface([ellipse_loop], tag=2)
occ.synchronize()
# --- Physical groups ---
gmsh.model.addPhysicalGroup(1, [4], tag=1) # left edge
gmsh.model.addPhysicalGroup(1, [1], tag=2) # top edge
gmsh.model.addPhysicalGroup(1, [2], tag=3) # right edge
gmsh.model.addPhysicalGroup(1, [3], tag=4) # bottom edge
gmsh.model.addPhysicalGroup(1, [8, 5, 6, 7], tag=5) # ellipse boundary
gmsh.model.addPhysicalGroup(2, [s_matrix], tag=1) # matrix
gmsh.model.addPhysicalGroup(2, [s_incl], tag=2) # inclusion
# Optional: add names
gmsh.model.setPhysicalName(2, 1, "Matrix")
gmsh.model.setPhysicalName(2, 2, "Inclusion")
gmsh.model.setPhysicalName(1, 5, "EllipseBoundary")
# --- Mesh ---
gmsh.option.setNumber("Mesh.Algorithm", 6)
gmsh.model.mesh.generate(2)
# gmsh.write("square_with_circle.msh")
model_rank = 0
gdim = 2
mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, model_rank, gdim)
# # gmsh.fltk.run() # Uncomment to view GUI
gmsh.finalize()
################# Meshing complete ##################################
################# Projection starts ###################################
res = "-"
degree = 1
V = dolfinx.fem.functionspace(mesh, ("CG", degree, (mesh.geometry.dim,)))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
print(f" Num dofs global {V.dofmap.index_map.size_global * V.dofmap.index_map_bs}")
# define integration measures
dS = Measure("dS", domain=mesh, subdomain_data=ft, subdomain_id=5)
# form is defined on the boundary => matrix wil have many zero diagonal entries
a_bnd = ufl.inner(u(res), v(res)) * dS
form = fem.form(a_bnd)
# Find all dofs associated to the facets we are integrating over
interface_facets = ft.find(5)
boundary_blocks = dolfinx.fem.locate_dofs_topological(
V, V.mesh.topology.dim - 1, interface_facets
)
# Set diagonal for everything not associated with a facet in sparsity pattern
imap = V.dofmap.index_map
all_blocks = np.arange(imap.size_local + imap.num_ghosts, dtype=np.int32)
zero_blocks = np.delete(all_blocks, boundary_blocks)
start = time.perf_counter()
sp = fem.create_sparsity_pattern(form)
sp.insert_diagonal(zero_blocks)
sp.finalize()
end = time.perf_counter()
print(f"Create sparsity pattern {end - start}")
start = time.perf_counter()
A = cpp.la.petsc.create_matrix(mesh.comm, sp)
A.setOption(A.Option.NEW_NONZERO_ALLOCATION_ERR, 1)
fem.petsc.assemble_matrix(A, form)
A.assemble(A.AssemblyType.FLUSH)
end = time.perf_counter()
print(f"Assemble boundary {end - start}")
def bc_diag(A):
bc = dolfinx.fem.dirichletbc(
dolfinx.fem.Constant(mesh, np.zeros(u.ufl_shape, dtype=PETSc.ScalarType)),
zero_blocks,
V,
)
start = time.perf_counter()
dolfinx.cpp.fem.petsc.insert_diagonal(A, V._cpp_object, [bc._cpp_object], 1.0)
end = time.perf_counter()
print(f"BC insert diag {end - start}s")
return bc
bc = bc_diag(A)
start = time.perf_counter()
A.assemble()
end = time.perf_counter()
print(f"Matrix communication {end - start}")
#
print(f"min(diag) = {A.getDiagonal().array.min()}")
print(f"number of diagonal zeros = {sum(A.getDiagonal().array == 0)}")
print(f"norm(A) = {A.norm()}")
normalVec = ufl.FacetNormal(mesh)
b = dolfinx.fem.form(ufl.inner(normalVec(res), v(res)) * dS)
b_vec = dolfinx.fem.petsc.assemble_vector(b)
dolfinx.fem.petsc.apply_lifting(b_vec, [form], [[bc]])
b_vec.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
dolfinx.fem.petsc.set_bc(b_vec, [bc])
ksp = PETSc.KSP().create(mesh.comm)
ksp.setOperators(A)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.setErrorIfNotConverged(True)
uh = dolfinx.fem.Function(V)
start = time.perf_counter()
ksp.solve(b_vec, uh.x.petsc_vec)
ksp.destroy()
b_vec.destroy()
end = time.perf_counter()
uh.x.scatter_forward()
print(f"CG solve {end - start}s")
with dolfinx.io.VTXWriter(mesh.comm, "normal_ellipse.bp", [uh], engine="BP4") as writer:
writer.write(0.0)