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!

9

I found that when describing the jump displacement of the interface (involving the problem of two subdomeshes): jump = u1(-)-u2(+), and its normal vector is expressed in the form of facetnormal(domain)(+). Since you said that there is no positive or negative distinction for the normal vectors of the internal lines in the sub-grid, then in which direction does the positive direction of the normal vector here point?
just like this code: Cohesive zone modeling restricted to an interface — Computational Mechanics Numerical Tours with FEniCSx

10

It clearly states in that demo what is + and what is -:

Moreover, we also consistently switch the orientation of the facets so that cells in subdomain 1 correspond to the "+" side of the interface and cells in subdomain 2 to the "-" side.

12

However, I came across a really intriguing phenomenon. In this demo, the jump displacement is described as u1(-) - u2(+) , and the facet normal is represented as ufl.Facetnormal(domain)(+) . Nevertheless, when it is changed to u1(+) - u2(-) with the facet normal being ufl.Facetnormal(domain)(-) , the outcome is completely different. I’m quite puzzled about how such a discrepancy occurs.That’s why i am stuck in the direction of the facetnormal

13

Note that, there is only one side for each displacement subdomain. For avoiding any issue with facet integration, the cells on the opposite side of a facet are mapped to the cell on the subdomain. Hence, when defining the jump u2("-")=u2("+") and u1("+")=u1("-"). So it is understandable that when changing the facet normal you change the complete orientation.

14

However, for the case of u1(+) - u2(-), the facet normal (-) corresponds to the outer normal vector of u2. And for u1(-) - u2(+), the facet normal (+) also corresponds to the outer normal vector of u2. According to your meaning, I understand that u1(+) - u2(-) is equal to u1(-) - u2(+). But when I change the positive and negative signs of u1 and u2, I also change the element to which the facet normal belongs. As a result, the facet normal always remains related to u2 and does not change. Unless we have specified the positive and negative signs of u1 and u2 from the very beginning, is that right?