I get the wrong one. Note that dx, dx_shape, dx_square have been tested carefully by integrating a function on them, and they are correct.
I set the f_shape and f_square profiles with
def interpolate_dg(f, g, sf=None, region_id=None):
Q = f.function_space()
element = Q.ufl_element()
if (element.family() != 'Discontinuous Lagrange'):
# the method has been called on a field defined on a continuous function space -> throw an error and exit
print(f'{col.Fore.RED}Error: interpolate_dg has been called on a field with a continuous function space!! Stopping now.{col.Fore.RESET}')
sys.exit(1)
if (element.value_shape() != g.value_shape()):
# the value shape of Q and that of g differ -> check whether this is due to a 'convention' issue where scalars have been given a value shape of (1,) vs. ()
if ((((element.value_shape() == ()) and (g.value_shape() == (1,))) or ((element.value_shape() == (1,)) and (g.value_shape() == ()))) == False):
# the discrepancy was not due to a convention issue -> throw an error and exit
print(f'{col.Fore.RED}Error: value shapes are different!!\n\telement value shape = {element.value_shape()}\n\tg value shape= {g.value_shape()}\nStopping now.{col.Fore.RESET}')
sys.exit(1)
value_size = int(np.prod(element.value_shape())) if element.value_shape() else 1
mesh = Q.mesh()
'''
dof_coordinates stores the coordinates of the points where DOFs sit. Because the field 'f' defined on each DOF has value_size components, dof_coordinates is composed of blocks, where each block has 'value_size' entries, and blocks are all identical
For example, dof_coordinates is of the form ->
row 0: [x0, y0] ← this corresponds to f[0] at DOF point 0
row 1: [x0, y0] ← this corresponds to f[1] at DOF point 0
...
row value_size [x1, y1] ← this corresponds to f[0] at DOF point 1
row value_size+1 [x1, y1] ← this corresponds to f[1] at DOF point 1
...
'''
dof_coordinates = Q.tabulate_dof_coordinates()
'''
f_values contains the value of 'f' on the DOFs, and it has the same structure as 'dof_coordinates'
f_values =
entry 0: f[0] at DOF point 0
entry 1: f[1] at DOF point 0
...
entry value_size f[0] at DOF point 1
entry value_size + 1 f[1] at DOF point 1
...
'''
f_values = f.vector().get_local() # get a copy of field values
for cell in cells(mesh):
# run on all mesh cells
if region_id != None:
#region_id has been given when calling this method -> evaluate sf on the cell to obtain the tag of cell `cell`
# compute 'sf' on the cell; under consideration
cell_tag = sf[cell]
'''
cell_dofs contains the IDs of the DOFs that are contained into 'cell', it has the structure
[
id_f_0_on_DOF_0,
id_f_0_on_DOF_1,
...,
id_f_0_on_DOF_{n_nodes-1},
id_f_1_on_DOF_0,
id_f_1_on_DOF_1,
...,
id_f_1_on_DOF_{n_nodes-1},
...
]
where the pattern is repeated value_size times, i.e., one for each component of 'f', and n_nodes = [number of DOFs in the cell] / [value_size]. In other words
cell_dofs[j * n_nodes + i] = [index in f.values().get_local() corresponding to the j-th component of the field 'f' sitting on ith DOF in the cell 'cell']
'''
cell_dofs = Q.dofmap().cell_dofs(cell.index())
n_nodes = len(cell_dofs) // value_size
'''
remove the redundancy in cell_dofs and store the result in
cell_dofs_unique =
[
id_DOF_0,
id_DOF_1,
...,
id_DOF_{n_nodes-1}
]
'''
cell_dofs_unique = cell_dofs[:n_nodes]
if (region_id == None) or (cell_tag == region_id):
# if 'cell_tag' == 'id', then the cell under consideration belongs to the surface tagged with 'id' -> set the DOFs of 'f' relative to this cell according to 'g'
for i in range(len(cell_dofs_unique)):
# run over physical DOFs contained to 'cell' and print out the value of 'f' by specifying that those DOFs belong to region tagged with 'cell_tag' in a separate column of the csv output file. Note that, because the space of 'f' is discontinuous, here DOFs in 'cell' may belong to different mesh regions, and thus have different tags
# pad 'x' to three dimensions
dof_coordinate = dof_coordinates[cell_dofs_unique[i]]
for j in range(value_size):
f_values[cell_dofs[j * n_nodes + i]] = np.atleast_1d(g(dof_coordinate[:2]))[j]
f.vector().set_local(f_values)
f.vector().apply("insert")
It’s very hard to give guidance with an example that only consist of snippets, and no reproducible code.
Could you try to reproduce this on a unit square where you mark a subdomain (for instance a square inside the unit square), and do the same thing? Then it should be possible to make a small reproducible example.
For instance, have you checked that interpolate_dg gives the expected functions when used in
i.e is fsp.f what you expect?
furthermore, you have not shown us how you define dx_shape and dx_square
from fenics import *
import ufl
import gmsh
import meshio
import numpy as np
# =========================================================
# 1. MESH GENERATION (Gmsh, with embedded interface)
# =========================================================
gmsh.initialize()
gmsh.model.add("model")
L, H = 1.0, 1.0
resolution = 0.05
shape_id = 1 # inner square (shape)
square_id = 2 # outer region
# outer boundary
p_out = [gmsh.model.geo.addPoint(x, y, 0)
for x, y in [(0,0),(L,0),(L,H),(0,H)]]
l_out = [gmsh.model.geo.addLine(p_out[i], p_out[(i+1)%4]) for i in range(4)]
outer_loop = gmsh.model.geo.addCurveLoop(l_out)
# inner shape (square)
shape_pts = [[0.3,0.3],[0.7,0.3],[0.7,0.7],[0.3,0.7]]
p_in = [gmsh.model.geo.addPoint(x, y, 0) for x, y in shape_pts]
l_in = [gmsh.model.geo.addLine(p_in[i], p_in[(i+1)%4]) for i in range(4)]
inner_loop = gmsh.model.geo.addCurveLoop(l_in)
surf_square = gmsh.model.geo.addPlaneSurface([outer_loop, inner_loop])
surf_shape = gmsh.model.geo.addPlaneSurface([inner_loop])
gmsh.model.geo.synchronize()
gmsh.model.mesh.embed(1, l_in, 2, surf_square) # embed interface in outer region
gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(2, [surf_shape], shape_id)
gmsh.model.addPhysicalGroup(2, [surf_square], square_id)
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", resolution)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", resolution)
gmsh.model.mesh.generate(2)
gmsh.write("mesh.msh")
gmsh.finalize()
# =========================================================
# 2. CONVERT TO FENICS (triangles only)
# =========================================================
def create_mesh(mesh_io, cell_type):
cells = mesh_io.get_cells_type(cell_type)
cell_data = mesh_io.get_cell_data("gmsh:physical", cell_type)
return meshio.Mesh(
points=mesh_io.points[:, :2],
cells=[(cell_type, cells)],
cell_data={"name_to_read": [cell_data]},
)
msh_io = meshio.read("mesh.msh")
meshio.write("mesh.xdmf", create_mesh(msh_io, "triangle"))
# =========================================================
# 3. LOAD MESH AND CELL FUNCTION
# =========================================================
mesh = Mesh()
with XDMFFile("mesh.xdmf") as fh:
fh.read(mesh)
sf = MeshFunction("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as fh:
fh.read(sf, "name_to_read")
# =========================================================
# 4. FACET CLASSIFICATION
# =========================================================
# We need three sets of interior facets:
# dS_shape : interface between shape and square subdomains
# dS_I_shape : interior facets within the shape subdomain
# dS_I_square : interior facets within the square subdomain
INTERFACE_ID = 3
INT_SHAPE_ID = 4
INT_SQUARE_ID = 5
mesh.init(mesh.topology().dim()-1, mesh.topology().dim()) # build facet-cell connectivity
mf = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
for facet in facets(mesh):
adj = list(cells(facet))
if len(adj) == 2:
c0, c1 = sf[adj[0]], sf[adj[1]]
if c0 != c1:
mf[facet] = INTERFACE_ID
elif c0 == shape_id:
mf[facet] = INT_SHAPE_ID
else:
mf[facet] = INT_SQUARE_ID
dx = Measure("dx", domain=mesh)
dx_shape = Measure("dx", domain=mesh, subdomain_data=sf, subdomain_id=shape_id)
dx_square = Measure("dx", domain=mesh, subdomain_data=sf, subdomain_id=square_id)
ds_outer = Measure("ds", domain=mesh) # outer boundary
dS_shape = Measure("dS", domain=mesh, subdomain_data=mf, # shape/square interface
subdomain_id=INTERFACE_ID)
dS_I_shape = Measure("dS", domain=mesh, subdomain_data=mf, # interior of shape
subdomain_id=INT_SHAPE_ID)
dS_I_square = Measure("dS", domain=mesh, subdomain_data=mf, # interior of square
subdomain_id=INT_SQUARE_ID)
# =========================================================
# 5. PLUS/MINUS DETERMINATION
# =========================================================
# Determine which UFL restriction side ('+' or '-') corresponds to each
# subdomain on the interface facets.
# In legacy FEniCS the '+' side is the cell that appears first in the
# facet-to-cell connectivity, which follows the global cell-index ordering.
# NOTE: this yields a *globally consistent* label only when all interface
# facets have the same subdomain on the same side – verify with the print.
def plus_minus(mesh, sf, id_0, id_1):
"""
Returns (label_id_0, label_id_1) ∈ {'+', '-'} by inspecting the first
interface facet. Assumes consistent orientation across all interface facets.
"""
mesh.init(mesh.topology().dim()-1, mesh.topology().dim())
for facet in facets(mesh):
adj = list(cells(facet))
if len(adj) == 2:
c0, c1 = sf[adj[0]], sf[adj[1]]
if {c0, c1} == {id_0, id_1}:
# FEniCS convention: the cell with lower global index → '+' side
plus_cell = adj[0] if adj[0].index() < adj[1].index() else adj[1]
if sf[plus_cell] == id_0:
return '+', '-'
else:
return '-', '+'
return '+', '-' # fallback
shape_label, square_label = plus_minus(mesh, sf, shape_id, square_id)
print(f'label_shape = {shape_label}')
print(f'label_square = {square_label}')
# =========================================================
# 6. FUNCTION SPACE
# =========================================================
Q = FunctionSpace(mesh, "DG", 2)
u = Function(Q)
nu_u = TestFunction(Q)
f = Function(Q)
u_exact = Function(Q)
# =========================================================
# 7. INTERPOLATION (subdomain-wise)
# =========================================================
def interpolate_dg(func, expr_func, sf, region_id=None):
coords = func.function_space().tabulate_dof_coordinates()
vals = func.vector().get_local()
for cell in cells(mesh):
if region_id is None or sf[cell] == region_id:
for dof in Q.dofmap().cell_dofs(cell.index()):
vals[dof] = expr_func(coords[dof])
func.vector().set_local(vals)
func.vector().apply("insert")
def f_shape(x): return 2.0 * (-1 + x[0]**2 + 10*x[1]**2)
def f_square(x): return 6.0
def u_sh(x): return 1 + x[0]**2 - 2*x[1]**2
def u_sq(x): return 1 + x[0]**2 + 2*x[1]**2
interpolate_dg(f, f_shape, sf, shape_id)
interpolate_dg(f, f_square, sf, square_id)
interpolate_dg(u_exact, u_sh, sf, shape_id)
interpolate_dg(u_exact, u_sq, sf, square_id)
# =========================================================
# 8. PARAMETERS AND GEOMETRIC QUANTITIES
# =========================================================
alpha = Constant(1e2)
n = FacetNormal(mesh)
h = CellDiameter(mesh) # NB: must be restricted ('.+') in dS integrals
# because CellEdgeVectors is discontinuous
i = ufl.indices(1)[0] # free UFL index (Einstein summation)
u.assign(Constant(1.0))
r_min = mesh.hmin()
# =========================================================
# 9. DG OPERATORS
# =========================================================
def average(expr):
"""Arithmetic average across an interior facet."""
return (expr('+') + expr('-')) / 2
def jump(f, n):
"""Vectorial jump of scalar f: f(+) n(+) + f(-) n(-)."""
return f('+')*n('+') + f('-')*n('-')
# =========================================================
# 10. VARIATIONAL FORM (exact transcription of original)
# =========================================================
# ---------- F_0 : bulk terms (from IBP) ----------
# shape region : nonlinear equation u * u_{,i} * (nu_u)_{,i} + f * nu_u
# square region : linear equation u_{,i} * (nu_u)_{,i} + f * nu_u
# Boundary/interface flux terms:
# outer boundary : - n_i u_{,i} nu_u
# interface, square side : - n_i u_{,i} nu_u (from linear BVP)
# interface, shape side : - n_i u u_{,i} nu_u (from nonlinear BVP)
def ufl_conditional_form(mesh, sf, form_a, form_b, tag_a, tag_b):
Q = FunctionSpace(mesh, 'DG', 0)
cell_tag = Function(Q)
cell_tag.vector()[:] = sf.array() # DG0 dofs are ordered by cell index
# build f as a pure UFL expression
result = conditional(
ufl.eq(cell_tag, tag_a),
form_a,
conditional(ufl.eq(cell_tag, tag_b),
form_b,
form_b * 0
)
)
return result
# option 1 (not working)
'''
F_0 = (
u * u.dx(i) * nu_u.dx(i) + f * nu_u
) * dx_shape + (
u.dx(i) * nu_u.dx(i) + f * nu_u
) * dx_square \
- n[i] * u.dx(i) * nu_u * ds_outer \
- n(square_label)[i] * (u(square_label)).dx(i) * nu_u(square_label) * dS_shape \
- n(shape_label)[i] * u(shape_label) * (u(shape_label)).dx(i) * nu_u(shape_label) * dS_shape
'''
# option 2 (working)
F_0 = ufl_conditional_form(mesh,
sf,
u * u.dx(i) * nu_u.dx(i) + f * nu_u,
u.dx(i) * nu_u.dx(i) + f * nu_u,
shape_id,
square_id
) * dx\
- n[i] * u.dx(i) * nu_u * ds_outer \
- n(square_label)[i] * (u(square_label)).dx(i) * nu_u(square_label) * dS_shape \
- n(shape_label)[i] * u(shape_label) * (u(shape_label)).dx(i) * nu_u(shape_label) * dS_shape
# ---------- F_I : interior DG (SIPG-style) ----------
# square interior: standard symmetric term -avg(grad u) · jump(nu_u, n)
# shape interior: nonlinear analog -avg(grad u) · jump(u nu_u, n)
# penalty on both: (alpha/h) jump(u,n) · jump(nu_u,n)
#
# NB: h must be restricted on dS integrals (CellEdgeVectors is discontinuous)
F_I = \
- average(u.dx(i)) * jump(nu_u, n)[i] * dS_I_square \
- average(u.dx(i)) * jump(u * nu_u, n)[i] * dS_I_shape \
+ alpha / r_min * (
jump(u, n)[i] * jump(nu_u, n)[i]
) * dS_I_square \
+ alpha / r_min * (
jump(u, n)[i] * jump(nu_u, n)[i]
) * dS_I_shape
# ---------- F_b : Nitsche boundary penalisation ----------
# outer boundary : (alpha/h) (u - u_exact) nu_u ds
# interface, square side : (alpha/h) (u - u_exact) nu_u dS
# interface, shape side : (alpha/h) (u - u_exact) nu_u dS
F_b = alpha / h * (u - u_exact) * nu_u * ds_outer \
+ alpha / r_min * (
(u(square_label) - u_exact(square_label)) * nu_u(square_label)
) * dS_shape \
+ alpha / r_min * (
(u(shape_label) - u_exact(shape_label)) * nu_u(shape_label)
) * dS_shape
F = F_0 + F_I + F_b
# =========================================================
# 11. SOLVE
# =========================================================
solve(F == 0, u)
print(f'error_shape = {assemble((u - u_exact)**2 * dx_shape)}')
print(f'error_square = {assemble((u - u_exact)**2 * dx_square)}')
The mesh is a square with a circle in it, where the circle is meshed inside. \Omega_\text{square} is the region between the circle and the square boundary, and \Omega_\text{O} the region inside the circle.
I solve \partial_i (u \,\partial_i u ) = f_\text{O} \text{ in } \Omega_\text{O}, \partial_i \partial_i u = f_\text{sq} \text{ in } \Omega_\text{sq}.
The code can run with two options:
Option 1: the two variational problems in \Omega_\text{O} and \Omega_\text{sq} are defind by using dx_shapedx_square, where each dx is multiplied by the respective functional.
Option 2: the two variational problem are implemented by using the global dx and switching between \Omega_\text{O} and \Omega_\text{sq} with ufl_conditional_form
Please uncomment / comment the code block relative to each option to run the code with that option.
The solution of Option 2 agrees with the exact one to within machine precision as it should (the exact solution is a second order polynomial and we use second-order spaces)
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.073e+03 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 3.098e-02 (tol = 1.000e-10) r (rel) = 2.886e-05 (tol = 1.000e-09)
Newton iteration 2: r (abs) = 2.035e-03 (tol = 1.000e-10) r (rel) = 1.895e-06 (tol = 1.000e-09)
Newton iteration 3: r (abs) = 8.184e-06 (tol = 1.000e-10) r (rel) = 7.625e-09 (tol = 1.000e-09)
Newton iteration 4: r (abs) = 1.631e-10 (tol = 1.000e-10) r (rel) = 1.520e-13 (tol = 1.000e-09)
Newton solver finished in 4 iterations and 4 linear solver iterations.
error_shape = 4.639144890041902e-22
error_square = 1.906535228505125e-28
The solution with Option 1 does not agree to within machine precision
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.073e+03 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 6.052e-01 (tol = 1.000e-10) r (rel) = 5.639e-04 (tol = 1.000e-09)
Newton iteration 2: r (abs) = 3.228e-03 (tol = 1.000e-10) r (rel) = 3.007e-06 (tol = 1.000e-09)
Newton iteration 3: r (abs) = 8.185e-06 (tol = 1.000e-10) r (rel) = 7.625e-09 (tol = 1.000e-09)
Newton iteration 4: r (abs) = 1.631e-10 (tol = 1.000e-10) r (rel) = 1.520e-13 (tol = 1.000e-09)
Newton solver finished in 4 iterations and 4 linear solver iterations.
error_shape = 4.639144890041902e-22
error_square = 6.106107010627693e-05
For more complicated problems, the solution of Option 1 is plain wrong. I tested thoroughly interpolate_dg.
A better way to restrict this is to follow the suggestion of @MiroK from:
which ensures a consistent orientation. I.e. mark all cells on the side you would like to have a “+” restriction on with a larger value than on the other side of the interface.
Instead of defining two different measures, I would suggest using a single measure, like
and then solve the problem twice to see if I can spot the issue:
a = assemble(F)
a2 = assemble(F2)
np.testing.assert_allclose(a.get_local(), a2.get_local())
J = derivative(F, u)
A = assemble(J)
A2 = assemble(derivative(F2, u))
np.testing.assert_allclose(A.array(), A2.array(), atol=1e-13)
solve(F2 == 0, u, solver_parameters={"newton_solver": {"linear_solver": "mumps"}})
d = u.copy(deepcopy=True)
u.vector()[:] = 1.0
solve(F == 0, u, solver_parameters={"newton_solver": {"linear_solver": "mumps"}})
and then I get that the jacobians are the same at the first iteration, but they seem to deviate afterwards, as I get the two solver outputs:
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.073e+03 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 3.098e-02 (tol = 1.000e-10) r (rel) = 2.886e-05 (tol = 1.000e-09)
Newton iteration 2: r (abs) = 2.035e-03 (tol = 1.000e-10) r (rel) = 1.895e-06 (tol = 1.000e-09)
Newton iteration 3: r (abs) = 8.184e-06 (tol = 1.000e-10) r (rel) = 7.625e-09 (tol = 1.000e-09)
Newton iteration 4: r (abs) = 1.631e-10 (tol = 1.000e-10) r (rel) = 1.520e-13 (tol = 1.000e-09)
Newton solver finished in 4 iterations and 4 linear solver iterations.
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 1.073e+03 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-09)
Newton iteration 1: r (abs) = 6.052e-01 (tol = 1.000e-10) r (rel) = 5.639e-04 (tol = 1.000e-09)
Newton iteration 2: r (abs) = 3.228e-03 (tol = 1.000e-10) r (rel) = 3.007e-06 (tol = 1.000e-09)
Newton iteration 3: r (abs) = 8.185e-06 (tol = 1.000e-10) r (rel) = 7.625e-09 (tol = 1.000e-09)
Newton iteration 4: r (abs) = 1.631e-10 (tol = 1.000e-10) r (rel) = 1.520e-13 (tol = 1.000e-09)
Newton solver finished in 4 iterations and 4 linear solver iterations.
where we see that the initial residual is the same.
However, if you lessen the restriction on the jacobians, you will get that:
Traceback (most recent call last):
File "/root/shared/mwe_mekong.py", line 286, in <module>
np.testing.assert_allclose(A.array(), A2.array())
File "/usr/lib/python3/dist-packages/numpy/testing/_private/utils.py", line 1530, in assert_allclose
assert_array_compare(compare, actual, desired, err_msg=str(err_msg),
File "/usr/lib/python3/dist-packages/numpy/testing/_private/utils.py", line 844, in assert_array_compare
raise AssertionError(msg)
AssertionError:
Not equal to tolerance rtol=1e-07, atol=0
Mismatched elements: 2753 / 34152336 (0.00806%)
Max absolute difference: 1.52100554e-14
Max relative difference: 82.
x: array([[37.001268, -4.704807, -4.403054, ..., 0. , 0. ,
0. ],
[-4.704807, 35.352108, -3.959018, ..., 0. , 0. ,...
y: array([[37.001268, -4.704807, -4.403054, ..., 0. , 0. ,
0. ],
[-4.704807, 35.352108, -3.959018, ..., 0. , 0. ,...
which is floating error accumulated in the conditional logic.
However, I think there might be a bug in FFC here, and I have no further time to play around with it. As I cannot reproduce the descrepancy in DOLFINx, this is likely something in legacy FFC:
from mpi4py import MPI
from dolfinx import default_scalar_type
from dolfinx.fem import functionspace, Function
from ufl import Measure, FacetNormal, CellDiameter, conditional, TestFunction, indices, eq, derivative
from dolfinx.io.gmsh import read_from_msh
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.fem import assemble_scalar, IntegralType, compute_integration_domains, Constant
from dolfinx.mesh import compute_incident_entities, meshtags, GhostMode
from scifem import compute_interface_data
import gmsh
import numpy as np
from dolfinx.cpp.mesh import create_cell_partitioner
# =========================================================
# 1. MESH GENERATION (Gmsh, with embedded interface)
# =========================================================
gmsh.initialize()
gmsh.model.add("model")
L, H = 1.0, 1.0
resolution = 0.05
shape_id = 1 # inner square (shape)
square_id = 2 # outer region
# outer boundary
p_out = [gmsh.model.geo.addPoint(x, y, 0)
for x, y in [(0,0),(L,0),(L,H),(0,H)]]
l_out = [gmsh.model.geo.addLine(p_out[i], p_out[(i+1)%4]) for i in range(4)]
outer_loop = gmsh.model.geo.addCurveLoop(l_out)
# inner shape (square)
shape_pts = [[0.3,0.3],[0.7,0.3],[0.7,0.7],[0.3,0.7]]
p_in = [gmsh.model.geo.addPoint(x, y, 0) for x, y in shape_pts]
l_in = [gmsh.model.geo.addLine(p_in[i], p_in[(i+1)%4]) for i in range(4)]
inner_loop = gmsh.model.geo.addCurveLoop(l_in)
surf_square = gmsh.model.geo.addPlaneSurface([outer_loop, inner_loop])
surf_shape = gmsh.model.geo.addPlaneSurface([inner_loop])
gmsh.model.geo.synchronize()
gmsh.model.mesh.embed(1, l_in, 2, surf_square) # embed interface in outer region
gmsh.model.geo.synchronize()
gmsh.model.addPhysicalGroup(2, [surf_shape], shape_id)
gmsh.model.addPhysicalGroup(2, [surf_square], square_id)
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", resolution)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", resolution)
gmsh.model.mesh.generate(2)
gmsh.write("mesh.msh")
gmsh.finalize()
gmode = GhostMode.shared_facet
partitioner = create_cell_partitioner(gmode, max_facet_to_cell_links=2)
mesh_data = read_from_msh("mesh.msh",MPI.COMM_WORLD, 0, gdim=2, partitioner=partitioner)
mesh = mesh_data.mesh
sf = mesh_data.cell_tags
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
num_facets_on_proc = mesh.topology.index_map(mesh.topology.dim-1).size_local + mesh.topology.index_map(mesh.topology.dim-1).num_ghosts
INT_SHAPE_ID = 1
INT_SQUARE_ID = 2
INT_INTERFACE_ID = INT_SHAPE_ID + INT_SQUARE_ID
facet_markers = np.zeros(num_facets_on_proc, dtype=np.int32)
facet_markers[compute_incident_entities(mesh.topology, sf.find(square_id), mesh.topology.dim, mesh.topology.dim-1)] += INT_SQUARE_ID
facet_markers[compute_incident_entities(mesh.topology, sf.find(shape_id), mesh.topology.dim, mesh.topology.dim-1)] += INT_SHAPE_ID
ft = meshtags(mesh, mesh.topology.dim-1, np.arange(num_facets_on_proc, dtype=np.int32), facet_markers)
interface_entities = compute_interface_data(sf, ft.find(INT_INTERFACE_ID))
integral_data = [(INT_INTERFACE_ID, interface_entities.flatten())]
integral_data +=[(INT_SHAPE_ID, compute_integration_domains(
IntegralType.interior_facet,
mesh.topology,
ft.find(shape_id),))]
integral_data +=[(INT_SQUARE_ID, compute_integration_domains(
IntegralType.interior_facet,
mesh.topology,
ft.find(square_id),))]
if shape_id > square_id:
shape_label = "-"
square_label = "+"
else:
shape_label = "+"
square_label = "-"
metadata = {"quadrature_degree": 15, "quadrature_scheme": "default"}
dx = Measure("dx", domain=mesh, subdomain_data=sf, metadata=metadata)
dx_shape = dx(shape_id)
dx_square = dx(square_id)
dS = Measure("dS", domain=mesh, subdomain_data=integral_data, metadata=metadata) # all interior facets
ds_outer = Measure("ds", domain=mesh, metadata=metadata) # outer boundary
dS_shape = dS(INT_INTERFACE_ID)
dS_I_shape = dS(INT_SHAPE_ID)
dS_I_square = dS(INT_SQUARE_ID)
if shape_id > square_id:
shape_label = "+"
square_label = "-"
else:
shape_label = "-"
square_label = "+"
print(f'label_shape = {shape_label} {shape_id=}')
print(f'label_square = {square_label} {square_id=}')
# =========================================================
# 6. FUNCTION SPACE
# =========================================================
Q = functionspace(mesh, ("DG", 2))
u = Function(Q)
nu_u = TestFunction(Q)
f = Function(Q)
u_exact = Function(Q)
# =========================================================
# 7. INTERPOLATION (subdomain-wise)
# =========================================================
def f_shape(x): return 2.0 * (-1 + x[0]**2 + 10*x[1]**2)
def f_square(x): return np.full_like(x[0], 6.0)
def u_sh(x): return 1 + x[0]**2 - 2*x[1]**2
def u_sq(x): return 1 + x[0]**2 + 2*x[1]**2
f.interpolate(f_shape, cells0=sf.find(shape_id))
f.interpolate(f_square, cells0=sf.find(square_id))
f.x.scatter_forward() # ensure f is up to date on all processes for assembly
u_exact.interpolate(u_sh, cells0=sf.find(shape_id))
u_exact.interpolate(u_sq, cells0=sf.find(square_id))
u_exact.x.scatter_forward() # ensure u_exact is up to date on all processes for assembly
# =========================================================
# 8. PARAMETERS AND GEOMETRIC QUANTITIES
# =========================================================
alpha = Constant(mesh,default_scalar_type(100))
n = FacetNormal(mesh)
h = CellDiameter(mesh) # NB: must be restricted ('.+') in dS integrals
# because CellEdgeVectors is discontinuous
i = indices(1)[0] # free UFL index (Einstein summation)
u.x.array[:] = 1.0 # initial guess (for computing h on dS)
h_loc = mesh.h(mesh.topology.dim, np.arange(mesh.topology.index_map(mesh.topology.dim).size_local, dtype=np.int32))
r_min = mesh.comm.allreduce(h_loc.min(), op=MPI.MIN)
# =========================================================
# 9. DG OPERATORS
# =========================================================
def average(expr):
"""Arithmetic average across an interior facet."""
return (expr('+') + expr('-')) / 2
def jump(f, n):
"""Vectorial jump of scalar f: f(+) n(+) + f(-) n(-)."""
return f('+')*n('+') + f('-')*n('-')
# =========================================================
# 10. VARIATIONAL FORM (exact transcription of original)
# =========================================================
# ---------- F_0 : bulk terms (from IBP) ----------
# shape region : nonlinear equation u * u_{,i} * (nu_u)_{,i} + f * nu_u
# square region : linear equation u_{,i} * (nu_u)_{,i} + f * nu_u
# Boundary/interface flux terms:
# outer boundary : - n_i u_{,i} nu_u
# interface, square side : - n_i u_{,i} nu_u (from linear BVP)
# interface, shape side : - n_i u u_{,i} nu_u (from nonlinear BVP)
def ufl_conditional_form(mesh, sf, form_a, form_b, tag_a, tag_b):
Q = functionspace(mesh, ("DG", 0))
cell_tag = Function(Q)
cell_tag.x.array[:] = sf.values # DG0 dofs are ordered by cell index
# build f as a pure UFL expression
result = conditional(
eq(cell_tag, tag_a),
form_a,
conditional(eq(cell_tag, tag_b),
form_b,
0.0
)
)
return result
# option 1 (not working)
i = indices(1)[0] # free UFL index (Einstein summation)
F_0 = (
u * u.dx(i) * nu_u.dx(i) + f * nu_u
) * dx_shape + (
u.dx(i) * nu_u.dx(i) + f * nu_u
) * dx_square
# option 2 (working)
F_01 = ufl_conditional_form(mesh,
sf,
u * u.dx(i) * nu_u.dx(i) + f * nu_u,
u.dx(i) * nu_u.dx(i) + f * nu_u,
shape_id,
square_id
) * dx
F_other = - n[i] * u.dx(i) * nu_u * ds_outer \
- n(square_label)[i] * (u(square_label)).dx(i) * nu_u(square_label) * dS_shape \
- n(shape_label)[i] * u(shape_label) * (u(shape_label)).dx(i) * nu_u(shape_label) * dS_shape
# ---------- F_I : interior DG (SIPG-style) ----------
# square interior: standard symmetric term -avg(grad u) · jump(nu_u, n)
# shape interior: nonlinear analog -avg(grad u) · jump(u nu_u, n)
# penalty on both: (alpha/h) jump(u,n) · jump(nu_u,n)
#
# NB: h must be restricted on dS integrals (CellEdgeVectors is discontinuous)
F_I = \
- average(u.dx(i)) * jump(nu_u, n)[i] * dS_I_square \
- average(u.dx(i)) * jump(u * nu_u, n)[i] * dS_I_shape \
+ alpha / r_min * (
jump(u, n)[i] * jump(nu_u, n)[i]
) * dS_I_square \
+ alpha / r_min * (
jump(u, n)[i] * jump(nu_u, n)[i]
) * dS_I_shape
# ---------- F_b : Nitsche boundary penalisation ----------
# outer boundary : (alpha/h) (u - u_exact) nu_u ds
# interface, square side : (alpha/h) (u - u_exact) nu_u dS
# interface, shape side : (alpha/h) (u - u_exact) nu_u dS
F_b = alpha / h * (u - u_exact) * nu_u * ds_outer \
+ alpha / r_min * (
(u(square_label) - u_exact(square_label)) * nu_u(square_label)
) * dS_shape \
+ alpha / r_min * (
(u(shape_label) - u_exact(shape_label)) * nu_u(shape_label)
) * dS_shape
F = F_0 + F_other + F_I + F_b
F2 = F_01 + F_other+ F_I + F_b
# =========================================================
# 11. SOLVE
# =========================================================
np.random.seed(0)
u.x.array[:] = np.random.rand(u.x.array.size) # random initial guess
from dolfinx.fem import assemble_vector, form, assemble_matrix
a = assemble_vector(form(F))
a2 = assemble_vector(form(F2))
np.testing.assert_allclose(a.array, a2.array)
J = derivative(F, u)
A = assemble_matrix(form(J))
A2 = assemble_matrix(form(derivative(F2, u)))
np.testing.assert_allclose(A.data, A2.data)
u.x.array[:] = 1.0
opts = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps", "snes_linesearch_type": "none", "snes_rtol": 1e-10,
"snes_error_if_not_converged": True, "snes_monitor": None}
solver = NonlinearProblem(F, u, petsc_options_prefix="org_solver", petsc_options=opts)
solver.solve()
d_array = u.x.array.copy()
u.x.array[:] = 1.0
solver_2 = NonlinearProblem(F2, u, petsc_options_prefix="option_solver", petsc_options=opts)
solver_2.solve()
print(np.max(np.abs(d_array - u.x.array)))
print(f'error_shape = {mesh.comm.allreduce(assemble_scalar(form((u - u_exact)**2 * dx_shape)), op=MPI.SUM)}')
print(f'error_square = {mesh.comm.allreduce(assemble_scalar(form((u - u_exact)**2 * dx_square)), op=MPI.SUM)}')
yielding:
label_shape = - shape_id=1
label_square = + square_id=2
0 SNES Function norm 1.073369114053e+03
1 SNES Function norm 4.162191337893e-01
2 SNES Function norm 2.028425869492e-03
3 SNES Function norm 8.088894392143e-06
4 SNES Function norm 1.585554376335e-10
0 SNES Function norm 1.073369114053e+03
1 SNES Function norm 4.162191337893e-01
2 SNES Function norm 2.028425869492e-03
3 SNES Function norm 8.088894392143e-06
4 SNES Function norm 1.585554376335e-10
0.0
error_shape = 8.189560766580326e-09
error_square = 8.465276585475065e-08
Thank you. Given that I am not going to switch to Dolfinx anytime soon, may you please let me know whether you see anything wrong with using Option 2 ?
They should be equivalent, even if option two is a massive overkill in terms of engineering, which can cause quadrature estimates to explode if not controlled.