Wrong FacetNormal vector on internal boundaries

dolfinx
1

Hi,

I have a 2D axi-symmetric problem in which I cannot get the expected behavior of the FacetNormal vector. Below is a simple MWE in which I consider a semi-circle in a rectangular box. I want to integrate the normal vector to the circle on the circle’s circumference. I evaluate it using FacetNormal vector as well as by defining the surface normal by hand. The latter gives correct result that I can easily work out on paper too. But using the FacetNormal vector I get numbers which do not make sense at all. Below is a simple MWE that works it out:

# This script contains an MWE to show that normal vector to the surface is not correct

import dolfinx, gmsh, mpi4py, ufl
from dolfinx.io import gmshio

#%% create geometry

sph_R = 0.3 # radius of sphere
xc, yc = 0.0, 0.9 # location of sphere's center on vertical axis

gmsh.initialize()
gmsh.clear()
occ = gmsh.model.occ
gdim = 2
# displace the circle along y-axis
circ = occ.addDisk(xc, yc, 0, sph_R, sph_R)
cut_box = occ.add_rectangle(xc-sph_R, yc-sph_R, 0, sph_R, 2*sph_R)
occ.synchronize()
half, _ = occ.cut([(gdim, circ)], [(gdim, cut_box)])
box = occ.add_rectangle(0, yc-2*sph_R, 0, 2*sph_R, 4*sph_R)
final = occ.fragment([(gdim, box)], half)
occ.synchronize()

all_doms = gmsh.model.getEntities(gdim)
for j, dom in enumerate(all_doms):
    gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j+1)

all_edges = gmsh.model.getEntities(gdim-1)
for j, edge in enumerate(all_edges):
    gmsh.model.addPhysicalGroup(edge[0], [edge[1]], edge[1])

gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.refine()

model_rank = 0
mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, model_rank, gdim)
# gmsh.fltk.run()
gmsh.finalize()

#%% Evaluate integrals

r, z = ufl.SpatialCoordinate(mesh) # let's call horizontal axis as r and vertical z
n = ufl.FacetNormal(mesh) # numerical normal vector
n_an = ufl.as_vector((r/sph_R, (z-yc)/sph_R)) # analytical normal vector
dS = ufl.Measure('dS', domain=mesh, subdomain_data=ft)

sph_bnd_idx = (2, 3)

integrand = dolfinx.fem.form(n_an[0]* dS(sph_bnd_idx))  # outer boundaries
val1 = dolfinx.fem.assemble.assemble_scalar(integrand)

integrand = dolfinx.fem.form(n('-')[0]* dS(sph_bnd_idx))  # outer boundaries
val2 = dolfinx.fem.assemble.assemble_scalar(integrand)
print("Line integral of analytical nr = {}, Line integral of numerical nr = {}".format(val1, val2))

integrand = dolfinx.fem.form(n_an[1]* dS(sph_bnd_idx))  # outer boundaries
val3 = dolfinx.fem.assemble.assemble_scalar(integrand)

integrand = dolfinx.fem.form(n('-')[1]* dS(sph_bnd_idx))  # outer boundaries
val4 = dolfinx.fem.assemble.assemble_scalar(integrand)
print("Line integral of analytical nz = {}, Line integral of numerical nz = {}".format(val3, val4))

I also notice that the result with FacetNormal nz (i.e. n[1] along vertical axis) actually corresponds to the correct answer I should have gotten from n[0] and vice versa. I am unsure if it is a coincidence.

Thanks a lot for help!

2

This is due to the fact that interior facet integrals have no consistent ordering. The (“+” and “-”) is arbitrary for any facet.
For the example you have in mind, this can be circumvented with using a one sided facet integral where you send in the integration entities consistently from one side.

# Computing a onesided integral over an interior facet with consistent orientation
# Enabled by: https://github.com/FEniCS/dolfinx/pull/2269
# Copyright 2023 Jørgen S. Dokken
# SPDX: MIT

import dolfinx
import numpy as np
import ufl
from mpi4py import MPI

mesh = dolfinx.mesh.create_unit_square(
    MPI.COMM_WORLD, 10, 10, ghost_mode=dolfinx.mesh.GhostMode.shared_facet)

# Create connectivties required for defining integration entities
tdim = mesh.topology.dim
fdim = tdim - 1
mesh.topology.create_entities(fdim)
mesh.topology.create_connectivity(fdim, tdim)
mesh.topology.create_connectivity(tdim, fdim)


# Get number of cells on process
cell_map = mesh.topology.index_map(tdim)
num_cells = cell_map.size_local + cell_map.num_ghosts

# Create markers for each size of the interface
cell_values = np.ones(num_cells, dtype=np.int32)
cell_values[dolfinx.mesh.locate_entities(
    mesh, tdim, lambda x: x[0] <= 0.5+1e-13)] = 2
ct = dolfinx.mesh.meshtags(mesh, tdim, np.arange(
    num_cells, dtype=np.int32), cell_values)

facet_map = mesh.topology.index_map(fdim)
num_facets = facet_map.size_local + facet_map.num_ghosts

# Create facet markers
facet_values = np.ones(num_facets, dtype=np.int32)
facet_values[dolfinx.mesh.locate_entities(
    mesh, fdim, lambda x: np.isclose(x[0], 0.5))] = 2
ft = dolfinx.mesh.meshtags(mesh, fdim, np.arange(
    num_facets, dtype=np.int32), facet_values)

# Give a set of facets marked with a value (in this case 2), get a consistent orientation for an interior integral
facets_to_integrate = ft.find(2)

f_to_c = mesh.topology.connectivity(fdim, tdim)
c_to_f = mesh.topology.connectivity(tdim, fdim)
# Compute integration entities for a single facet of a cell.
# Each facet is represented as a tuple (cell_index, local_facet_index), where cell_index is local to process
# local_facet_index is the local indexing of a facet for a given cell
integration_entities = []
for i, facet in enumerate(facets_to_integrate):
    # Only loop over facets owned by the process to avoid duplicate integration
    if facet >= facet_map.size_local:
        continue
    # Find cells connected to facet
    cells = f_to_c.links(facet)
    # Get value of cells
    marked_cells = ct.values[cells]
    # Get the cell marked with 2
    correct_cell = np.flatnonzero(marked_cells == 2)

    assert len(correct_cell) == 1
    # Get local index of facet
    local_facets = c_to_f.links(cells[correct_cell[0]])
    local_index = np.flatnonzero(local_facets == facet)
    assert len(local_index) == 1

    # Append integration entities
    integration_entities.append(cells[correct_cell[0]])
    integration_entities.append(local_index[0])

# Create custom integration measure for one-sided integrals
ds = ufl.Measure("ds", domain=mesh, subdomain_data=[
                 (8, np.asarray(integration_entities, dtype=np.int32))])
n = ufl.FacetNormal(mesh)
x = ufl.SpatialCoordinate(mesh)
# Exact integral is [y/2**2]_0^1= 1/2
L = ufl.dot(ufl.as_vector((x[1], 0)), n)*ds(8)
L_compiled = dolfinx.fem.form(L)

print(
    f"Correct integral: {mesh.comm.allreduce(dolfinx.fem.assemble_scalar(L_compiled), op=MPI.SUM)}")


# Create reference implementation where we use a restricted two-sided integral with no notion of orientation
dS = ufl.Measure("dS", domain=mesh, subdomain_data=ft, subdomain_id=2)
n = ufl.FacetNormal(mesh)
x = ufl.SpatialCoordinate(mesh)
L2 = ufl.dot(ufl.as_vector((x[1], 0)), n("+"))*dS
L2_compiled = dolfinx.fem.form(L2)
print(
    f"Wrong integral: {mesh.comm.allreduce(dolfinx.fem.assemble_scalar(L2_compiled), op=MPI.SUM)}")

If you run this on 1,2,3 or 4 processes, you will observe that only the first approach gives a consistent and correct answer.

3 Likes
3

Thanks a lot for detailed response! I will need to adapt this to my actual problem. Like the normal vector, I guess we can work with any solution variable?

For now, could you please have a look again at this which is running into error on my end:

# Create custom integration measure for one-sided integrals
ds = ufl.Measure("ds", domain=mesh, subdomain_data=[
                 (8, np.asarray(integration_entities, dtype=np.int32))])
4

Yes, you can use anything that is in ufl.

What version of dolfinx are you running?

5

What version of dolfinx are you running?

Sorry I should have given more details. The error comes from the np.asarray function:

ValueError                                Traceback (most recent call last)
File ~/work/upwork/2023/magnetic_harvesting/fenics_models/dokken_integrate_example.py:77
     73     integration_entities.append(local_index)
     75 # Create custom integration measure for one-sided integrals
     76 ds = ufl.Measure("ds", domain=mesh, subdomain_data=[
---> 77                  (8, np.asarray(integration_entities, dtype=np.int32))])
     78 n = ufl.FacetNormal(mesh)
     79 x = ufl.SpatialCoordinate(mesh)

ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (20,) + inhomogeneous part.

I am using dolfinx 0.7.1 and numpy 1.26.2

6

I cannot reproduce this with the docker image: ghcr.io/fenics/dolfinx/dolfinx:v0.7.1
Make sure to post your full code.

7

Could you try changing the line:

to
integration_entities.append(local_index[0])
as I believe this is what causes the value error with modern numpy

8

It’s working now, thanks a lot for your help!