It is because the input to C++ requires integration_entities to be flattened (i.e. a flat array where each row is stacked after eachother).
1 Like
OK your example function works just fine, also on my mesh. However, as soon as I try to compute a tangential derivative, I get weird results, e.g. consider this code
import gmsh
from mpi4py import MPI
from dolfinx.io import gmshio
import meshio
import numpy as np
from ufl import Measure, grad, inner, FacetNormal, dot, SpatialCoordinate, as_vector
from dolfinx.fem import form, assemble_scalar, functionspace, Function
from dolfinx.mesh import meshtags
from dolfinx import plot
import pyvista as pv
def grad_t(u, n):
return grad(u) - dot(grad(u), n) * n
def get_mesh(size: float = 0.2, meshsize: float = 0.05, name: str = "dev"):
comm = MPI.COMM_WORLD
gmsh.initialize()
gmsh.model.add(name)
model = gmsh.model
geom = model.occ
gdim = 2
defeatured_domain_tag = geom.addDisk(0, 0, 0, 1, 1)
geom.synchronize()
disk_tag = geom.addDisk(0, 0, 0, size, size, defeatured_domain_tag + 1)
geom.synchronize()
out_dim_tags, *_ = geom.cut(
[(gdim, defeatured_domain_tag)],
[(gdim, disk_tag)],
removeTool=False,
removeObject=False,
)
_, exact_domain_tag = out_dim_tags[0]
geom.synchronize()
exact_domain_marker = model.addPhysicalGroup(gdim, [exact_domain_tag])
model.setPhysicalName(gdim, exact_domain_marker, "ExactDomain")
outer_boundary = model.getBoundary([(gdim, defeatured_domain_tag)], oriented=False)
outer_boundary_marker = model.addPhysicalGroup(
gdim - 1, [b[1] for b in outer_boundary]
)
model.setPhysicalName(gdim - 1, outer_boundary_marker, "OuterBoundary")
disk_marker = model.addPhysicalGroup(gdim, [disk_tag])
model.setPhysicalName(gdim, disk_marker, "Disk")
disk_boundary = model.getBoundary([(gdim, disk_tag)], oriented=False)
disk_boundary_marker = model.addPhysicalGroup(
gdim - 1, [b[1] for b in disk_boundary]
)
model.setPhysicalName(gdim - 1, disk_boundary_marker, "BoundaryDisk")
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", meshsize)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", meshsize)
gmsh.model.mesh.generate(dim=gdim)
gmsh.write(f"{name}.msh")
mesh = meshio.read(f"{name}.msh")
meshio.write(f"{name}.xdmf", mesh)
domain, cell_markers, facet_markers = gmshio.model_to_mesh(
gmsh.model, comm, rank=0, gdim=gdim
)
gmsh.finalize()
name2marker = {
"ExactDomain": exact_domain_marker,
"OuterBoundary": outer_boundary_marker,
"Disk": disk_marker,
"BoundaryDisk": disk_boundary_marker,
}
return domain, cell_markers, facet_markers, name2marker
if __name__ == "__main__":
size = 0.5
mesh, cell_tags, facet_tags, name2marker = get_mesh(size=size, meshsize=0.01)
tdim = mesh.topology.dim
fdim = tdim - 1
disk_boundary_marker = name2marker["BoundaryDisk"]
disk_marker = name2marker["Disk"]
outer_marker = name2marker["OuterBoundary"]
outer_facets = facet_tags.find(outer_marker)
outer_tags = meshtags(mesh, fdim, outer_facets, outer_marker)
disk_boundary_facets = facet_tags.find(disk_boundary_marker)
f2c = mesh.topology.connectivity(fdim, tdim)
c2f = mesh.topology.connectivity(tdim, fdim)
integration_entities = []
for i, facet in enumerate(disk_boundary_facets):
cells = f2c.links(facet)
marked_cells = cell_tags.values[cells]
correct_cell = np.flatnonzero(marked_cells == disk_marker)
assert len(correct_cell) == 1
local_facets = c2f.links(cells[correct_cell[0]])
local_index = np.flatnonzero(local_facets == facet)
assert len(local_index) == 1
integration_entities.append((cells[correct_cell[0]], local_index[0]))
ds = Measure("ds", domain=mesh, subdomain_data=[
(disk_boundary_marker, np.asarray(integration_entities, dtype=np.int32).flatten())
])
n = FacetNormal(mesh)
V = functionspace(mesh, ("CG", 1))
u = Function(V)
u.interpolate(lambda x: x[0]**2)
integral = inner(grad_t(u, n), grad_t(u, n)) * ds
norm = form(integral)
print(
f"Integral: {mesh.comm.allreduce(assemble_scalar(norm), op=MPI.SUM)}"
)
print(
f"Correct integral: {size**3 * np.pi}"
)
What would be the correct way to compute the tangential derivative of a function and then compute its L^2-norm?
What do you mean by weird results?
Could you also consider this on a simple, circular geometry with a standard measure?
Secondly, you are approximating x^2 with a piecewise linear function, which means that \nabla u is a piecewise constant, which would give you way lower accuracy than you would hope for. Try the same example with u being in a second order Lagrange space.
By weird results I mean that the integral of the tangential gradient on an interior boundary is not correctly computed, even if I explicitly construct an oriented measure as you have explained above. Here, by interior boundary, I mean the inner circle of my two circles.
When I do the same but only on the inner circle with the default measure on outer boundaries provided by FenicsX, the result is correct.
Here’s an MWE:
import gmsh
from mpi4py import MPI
from dolfinx.io import gmshio
import meshio
import numpy as np
from ufl import Measure, grad, inner, FacetNormal, dot
from dolfinx.fem import form, assemble_scalar, functionspace, Function
def grad_t(u, n):
return grad(u) - dot(grad(u), n) * n
def get_inner_mesh(size: float = 0.2, meshsize: float = 0.05, name: str = "dev_inner"):
comm = MPI.COMM_WORLD
gmsh.initialize()
gmsh.model.add(name)
model = gmsh.model
geom = model.occ
gdim = 2
domain_tag = geom.addDisk(0, 0, 0, size, size)
geom.synchronize()
exact_domain_marker = model.addPhysicalGroup(gdim, [domain_tag])
model.setPhysicalName(gdim, exact_domain_marker, "ExactDomain")
boundary = model.getBoundary([(gdim, domain_tag)], oriented=False)
boundary_marker = model.addPhysicalGroup(gdim - 1, [b[1] for b in boundary])
model.setPhysicalName(gdim - 1, boundary_marker, "OuterBoundary")
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", meshsize)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", meshsize)
gmsh.model.mesh.generate(dim=gdim)
gmsh.write(f"{name}.msh")
mesh = meshio.read(f"{name}.msh")
meshio.write(f"{name}.xdmf", mesh)
domain, cell_markers, facet_markers = gmshio.model_to_mesh(
gmsh.model, comm, rank=0, gdim=gdim
)
gmsh.finalize()
name2marker = {"ExactDomain": exact_domain_marker, "OuterBoundary": boundary_marker}
return domain, cell_markers, facet_markers, name2marker
def get_full_mesh(size: float = 0.2, meshsize: float = 0.05, name: str = "dev"):
comm = MPI.COMM_WORLD
gmsh.initialize()
gmsh.model.add(name)
model = gmsh.model
geom = model.occ
gdim = 2
defeatured_domain_tag = geom.addDisk(0, 0, 0, 1, 1)
geom.synchronize()
disk_tag = geom.addDisk(0, 0, 0, size, size, defeatured_domain_tag + 1)
geom.synchronize()
out_dim_tags, *_ = geom.cut(
[(gdim, defeatured_domain_tag)],
[(gdim, disk_tag)],
removeTool=False,
removeObject=False,
)
_, exact_domain_tag = out_dim_tags[0]
geom.synchronize()
exact_domain_marker = model.addPhysicalGroup(gdim, [exact_domain_tag])
model.setPhysicalName(gdim, exact_domain_marker, "ExactDomain")
outer_boundary = model.getBoundary([(gdim, defeatured_domain_tag)], oriented=False)
outer_boundary_marker = model.addPhysicalGroup(
gdim - 1, [b[1] for b in outer_boundary]
)
model.setPhysicalName(gdim - 1, outer_boundary_marker, "OuterBoundary")
disk_marker = model.addPhysicalGroup(gdim, [disk_tag])
model.setPhysicalName(gdim, disk_marker, "Disk")
disk_boundary = model.getBoundary([(gdim, disk_tag)], oriented=False)
disk_boundary_marker = model.addPhysicalGroup(
gdim - 1, [b[1] for b in disk_boundary]
)
model.setPhysicalName(gdim - 1, disk_boundary_marker, "BoundaryDisk")
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", meshsize)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", meshsize)
gmsh.model.mesh.generate(dim=gdim)
gmsh.write(f"{name}.msh")
mesh = meshio.read(f"{name}.msh")
meshio.write(f"{name}.xdmf", mesh)
domain, cell_markers, facet_markers = gmshio.model_to_mesh(
gmsh.model, comm, rank=0, gdim=gdim
)
gmsh.finalize()
name2marker = {
"ExactDomain": exact_domain_marker,
"OuterBoundary": outer_boundary_marker,
"Disk": disk_marker,
"BoundaryDisk": disk_boundary_marker,
}
return domain, cell_markers, facet_markers, name2marker
def get_oriented_measure(
mesh, cell_tags, tdim, fdim, int_boundary_marker, int_marker, int_boundary_facets
):
f2c = mesh.topology.connectivity(fdim, tdim)
c2f = mesh.topology.connectivity(tdim, fdim)
integration_entities = []
for i, facet in enumerate(int_boundary_facets):
cells = f2c.links(facet)
marked_cells = cell_tags.values[cells]
correct_cell = np.flatnonzero(marked_cells == int_marker)
assert len(correct_cell) == 1
local_facets = c2f.links(cells[correct_cell[0]])
local_index = np.flatnonzero(local_facets == facet)
assert len(local_index) == 1
integration_entities.append((cells[correct_cell[0]], local_index[0]))
ds = Measure(
"ds",
domain=mesh,
subdomain_data=[
(
int_boundary_marker,
np.asarray(integration_entities, dtype=np.int32).flatten(),
)
],
)
return ds
def run_full_mwe(size: float = 0.5, meshsize: float = 0.01, deg: int = 1):
print("Running MWE on the full disk with oriented measure...")
mesh, cell_tags, facet_tags, name2marker = get_full_mesh(
size=size, meshsize=meshsize
)
tdim = mesh.topology.dim
fdim = tdim - 1
disk_boundary_marker = name2marker["BoundaryDisk"]
disk_marker = name2marker["Disk"]
disk_boundary_facets = facet_tags.find(disk_boundary_marker)
ds = get_oriented_measure(
mesh,
cell_tags,
tdim,
fdim,
disk_boundary_marker,
disk_marker,
disk_boundary_facets,
)
n = FacetNormal(mesh)
V = functionspace(mesh, ("CG", 1))
u = Function(V)
u.interpolate(lambda x: x[0] ** 2)
integral = inner(grad_t(u, n), grad_t(u, n)) * ds
norm = form(integral)
result_title = "----- Result on full disk and interior boundary -----"
print(result_title)
print(f"Integral: {mesh.comm.allreduce(assemble_scalar(norm), op=MPI.SUM)}")
print(f"Correct integral: {size**3 * np.pi}")
print("-" * len(result_title))
def run_inner_mwe(size: float = 0.5, meshsize: float = 0.01, deg: int = 1):
print("Running MWE on the inner disk with default measure...")
mesh, *_ = get_inner_mesh(
size=size, meshsize=meshsize
)
ds = Measure("ds", domain=mesh)
n = FacetNormal(mesh)
V = functionspace(mesh, ("CG", 1))
u = Function(V)
u.interpolate(lambda x: x[0] ** 2)
integral = inner(grad_t(u, n), grad_t(u, n)) * ds
norm = form(integral)
result_title = "----- Result on inner disk -----"
print(result_title)
print(f"Integral: {mesh.comm.allreduce(assemble_scalar(norm), op=MPI.SUM)}")
print(f"Correct integral: {size**3 * np.pi}")
print("-" * len(result_title))
if __name__ == "__main__":
run_full_mwe(size=0.5, meshsize=0.005, deg=1)
run_inner_mwe(size=0.5, meshsize=0.005, deg=1)