Consistent side of interior boundaries with dolfinx

dolfinx
1

Hi

I want to create an interior discontinuity for Poisson’s equation as below

import dolfinx
import dolfinx.fem.petsc
from petsc4py import PETSc
import ufl
from mpi4py import MPI
import numpy as np

domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 8, 8)
V = dolfinx.fem.functionspace(domain, ("DG", 1))
tdim = domain.topology.dim
fdim = tdim - 1

def on_center_center(x):
    return np.isclose(x[0],0.5)&(x[1]>=.5)

def left_domain(x):
    return (x[0]<=.5)

domain.topology.create_connectivity(fdim, tdim)
domain.topology.create_connectivity(tdim, fdim)

facet_map = domain.topology.index_map(fdim)
num_facets_on_proc = facet_map.size_local + facet_map.num_ghosts
facets = np.arange(num_facets_on_proc, dtype=np.int32)
facet_values = np.ones(num_facets_on_proc, dtype=np.int32)
facet_values[dolfinx.mesh.locate_entities(domain, domain.topology.dim - 1, on_center_center)] = 2
ft = dolfinx.mesh.meshtags(domain, domain.topology.dim - 1, facets, facet_values)
dS = ufl.Measure("dS", domain=domain, subdomain_data=ft)

#Setting subdomain_data to define +/- side https://fenicsproject.discourse.group/t/integrating-over-an-interior-surface/247/4
cell_map = domain.topology.index_map(tdim)
num_cells_on_proc = cell_map.size_local + cell_map.num_ghosts
cells = np.arange(num_cells_on_proc, dtype=np.int32)
cell_values = 2*np.ones(num_cells_on_proc, dtype=np.int32)
cell_values[dolfinx.mesh.locate_entities(domain, tdim, left_domain)] = 1 # 3
st = dolfinx.mesh.meshtags(domain, tdim, cells, cell_values)
dx = ufl.Measure('dx', domain=domain,subdomain_data=st)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = dolfinx.fem.Constant(domain, dolfinx.default_scalar_type(0))
leftValue=dolfinx.fem.Constant(domain, dolfinx.default_scalar_type(-1))
rightValue=dolfinx.fem.Constant(domain, dolfinx.default_scalar_type(1))

n=ufl.FacetNormal(domain)
h = ufl.Circumradius(domain)
alpha = dolfinx.fem.Constant(domain, dolfinx.default_scalar_type(100))
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * dx   

a += -ufl.inner(ufl.jump(v,n),ufl.avg(ufl.grad(u)))*dS(1)
a += -ufl.inner(ufl.avg(ufl.grad(v)),ufl.jump(u,n))*dS(1)
a +=  alpha/ufl.avg(h)*ufl.inner(ufl.jump(v,n),ufl.jump(u,n)) * dS(1)

a += - ufl.inner(n('+'), ufl.grad(u('+'))) * v('+') * (dS(2))
a += - ufl.inner(n('+'), ufl.grad(v('+'))) * u('+') * (dS(2)) + alpha/h('+')  * ufl.inner(u('+'), v('+')) * (dS(2))

a += - ufl.inner(n('-'), ufl.grad(u('-'))) * v('-') * (dS(2))
a += - ufl.inner(n('-'), ufl.grad(v('-'))) * u('-') * (dS(2)) + alpha/h('-')  * ufl.inner(u('-'), v('-')) * (dS(2))

L = f * v * dx
L += - ufl.inner(n('+'), ufl.grad(v('+'))) * leftValue * dS(2) + alpha/h('+')  * ufl.inner(leftValue, v('+')) * dS(2)
L += - ufl.inner(n('-'), ufl.grad(v('-'))) * rightValue * dS(2) + alpha/h('-')  * ufl.inner(rightValue, v('-')) * dS(2)

A = dolfinx.fem.petsc.create_matrix(dolfinx.fem.form(a))
b = dolfinx.fem.petsc.create_vector(dolfinx.fem.form(L))
dolfinx.fem.petsc.assemble_vector(b, dolfinx.fem.form(L))
dolfinx.fem.petsc.assemble_matrix(A, dolfinx.fem.form(a))
A.assemble()
ksp=PETSc.KSP()
solver = ksp.create(domain.comm)
solver.setOperators(A)
uhd = dolfinx.fem.Function(V)
solver.solve(b, uhd.vector)

This works, but except for this simple case the (‘+’) and (‘-’) side is arbitrarily defined. I tried to apply this, but it does not seem to work with dolfinx.

How can consistently define the direction of “+”?

Thanks
Torben

2

You can follow the logic in: Wrong FacetNormal vector on internal boundaries - #2 by dokken
i.e. you pack your integration data by hand, and choose consistently that a cell marked with a value (say a will be the first in the integration entities array).

Please note that is probably better to use a dg-0 function to indicate right and left value.

Here is a reference implementation that does this handling consistently:

# Consistent interior integral
# Author: Jørgen S. Dokken
# SPDX-License-Identifier: MIT

from mpi4py import MPI
import dolfinx
import numpy as np
import packaging.version

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

# Divide mesh in two parts
def left(x):
    return x[0] < 0.5 + 1e-14

num_cells = mesh.topology.index_map(mesh.topology.dim).size_local+mesh.topology.index_map(mesh.topology.dim).num_ghosts
cells = np.arange(num_cells, dtype=np.int32)
values = np.full_like(cells, 3, dtype=np.int32) # Set all cells to 3
values[dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, left)] = 5 # Set left side to 5

# Create meshtag
ct = dolfinx.mesh.meshtags(mesh, mesh.topology.dim, cells, values)

# Find common facets
mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)
cell_left = ct.find(5)
cell_right = ct.find(3)
left_facets= dolfinx.mesh.compute_incident_entities(mesh.topology, cell_left, mesh.topology.dim, mesh.topology.dim-1)
right_facets= dolfinx.mesh.compute_incident_entities(mesh.topology, cell_right, mesh.topology.dim, mesh.topology.dim-1)
common_facets = np.intersect1d(left_facets, right_facets)


if packaging.version.parse(dolfinx.__version__) < packaging.version.parse("0.9"):

    integration_entities = []
    c_to_f = mesh.topology.connectivity(mesh.topology.dim, mesh.topology.dim - 1)
    f_to_c = mesh.topology.connectivity(mesh.topology.dim-1, mesh.topology.dim)

    for f in common_facets:
        cells = f_to_c.links(f)
        if len(cells) == 2:
            cell_values = ct.values[cells]
            if len(np.unique(cell_values)) < 2:
                # Only integrate over common boundary
                continue
            # We set all highest values first ("+") restriction
            argsort = np.argsort(cell_values)
            for cell in cells[argsort[::-1]]:
                facets = c_to_f.links(cell)
                pos = np.argwhere(facets == f)
                integration_entities.append(cell)
                integration_entities.append(pos[0, 0])
        else:
            pass
else:
    mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
    mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)
    integration_entities = dolfinx.fem.compute_integration_domains(dolfinx.fem.IntegralType.interior_facet, mesh.topology, common_facets,
                                                                   mesh.topology.dim-1)
    connected_cells = integration_entities.reshape(-1, 4)[:, [0, 2]]
    switch = ct.values[connected_cells[:,0]] < ct.values[connected_cells[:,1]]

    if len(switch)> 0 and  any(switch):
        tmp_entities = integration_entities.reshape(-1, 4)
        tmp_entities[switch] = tmp_entities[switch][:, [2, 3, 0, 1]]

        integration_entities = tmp_entities.flatten()

import ufl
dS_custom = ufl.Measure("dS", domain=mesh, subdomain_data=[(12, np.array(integration_entities, dtype=np.int32))])

# Reference implementation, which doesn't work
# facet_marker = dolfinx.mesh.meshtags(mesh, mesh.topology.dim-1, common_facets, np.full_like(common_facets, 12))
# dS_custom = ufl.Measure("dS", subdomain_data=facet_marker)


V = dolfinx.fem.functionspace(mesh, ("DG", 0, (2, )))
u = dolfinx.fem.Function(V)
# Set all values on right side to (2.5, 2.5)

if packaging.version.parse(dolfinx.__version__) < packaging.version.parse("0.9"):
    u.interpolate(lambda x: np.full((2, x.shape[1]), 2.5, dtype=np.float64), cells=ct.find(3))
else:
    u.interpolate(lambda x: np.full((2, x.shape[1]), 2.5, dtype=np.float64), cells0=ct.find(3))
u.x.scatter_forward()

# Set all values on left side to (0.3, 0)
def g(x):
    values = np.zeros((2, x.shape[1]), dtype=np.float64)
    values[0] = 0.3
    return values
if packaging.version.parse(dolfinx.__version__) < packaging.version.parse("0.9"):
    u.interpolate(g, cells=ct.find(5))
else:
    u.interpolate(g, cells0=ct.find(5))

u.x.scatter_forward()

# Evaluate integral from either side
n = ufl.FacetNormal(mesh)
f_left = dolfinx.fem.form(ufl.dot(u("+"),n("+")) * dS_custom(12))
left_val = mesh.comm.allreduce(dolfinx.fem.assemble_scalar(f_left), op=MPI.SUM)




f_right = dolfinx.fem.form(ufl.dot(u("-"),n("-")) * dS_custom(12))
right_val = mesh.comm.allreduce(dolfinx.fem.assemble_scalar(f_right), op=MPI.SUM)
print(left_val, right_val)
assert np.isclose(left_val, 0.3)
assert np.isclose(right_val, -2.5)
3

How do I define multiple integration domains using this method?

I have tried the following (and various other), but it does not seem to work.

dS_custom = ufl.Measure("dS", domain=domain,  subdomain_data=[(
    8, np.array(integration_entities, dtype=np.int32),
    10, np.array(integration_entities2, dtype=np.int32)
)])

And, is there a way to add “remaining” domains?

Thanks

4

Thanks the code below seems to work.

@dokken what do you mean DG-0 should work better? Instead of “Constant”?

What I’m actually trying to do is what in Comsol is called “Coil Geometry Analysis”, where they calculate the current distribution in a closed structure. Does anybody know how they do that?

import dolfinx
import dolfinx.fem.petsc
from petsc4py import PETSc
import ufl
from mpi4py import MPI
import numpy as np

domain = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 8, 8)
V = dolfinx.fem.functionspace(domain, ("DG", 1))
tdim = domain.topology.dim
fdim = tdim - 1

def right_domain(x):
    return (x[0]>=.75)

def inside_domain(x):
    return (x[0]<=.75)&(x[0]>=.5)&(x[1]>=.5)

def outside_domain(x):
    return (x[0]<=.75)&(x[0]>=.5)&(x[1]<=.5)

def left_domain(x):
    return (x[0]<=.5)

def internal_side_order(mesh,facets,cell_order):
    mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
    mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim - 1)
    integration_entities = dolfinx.fem.compute_integration_domains(dolfinx.fem.IntegralType.interior_facet, mesh.topology, facets,mesh.topology.dim-1)
    entities_tmp=integration_entities.reshape(-1, 4)
    switch=np.less(cell_order[entities_tmp[:,0]],cell_order[entities_tmp[:,2]])
    if len(switch)> 0 and  any(switch):
        entities_tmp[switch,:] = entities_tmp[switch,:][:, [2, 3, 0, 1]]
        integration_entities = entities_tmp.flatten()
    return(integration_entities)

num_cells = domain.topology.index_map(domain.topology.dim).size_local+domain.topology.index_map(domain.topology.dim).num_ghosts
cells = np.arange(num_cells, dtype=np.int32)
values = np.empty_like(cells, dtype=np.int32)

values[dolfinx.mesh.locate_entities(domain, domain.topology.dim, inside_domain)] = 2 # Set left side to 5
values[dolfinx.mesh.locate_entities(domain, domain.topology.dim, right_domain)] = 1 # Set left side to 5
values[dolfinx.mesh.locate_entities(domain, domain.topology.dim, outside_domain)] = 1 # Set left side to 5
values[dolfinx.mesh.locate_entities(domain, domain.topology.dim, left_domain)] = 3 # Set left side to 5

ct = dolfinx.mesh.meshtags(domain, domain.topology.dim, cells, values)

domain.topology.create_connectivity(fdim, tdim)
domain.topology.create_connectivity(tdim, fdim)

cell_inside = ct.find(1)
cell_outside = ct.find(2)
cell_left = ct.find(3)
inside_facets= dolfinx.mesh.compute_incident_entities(domain.topology, cell_inside, domain.topology.dim, domain.topology.dim-1)
outside_facets= dolfinx.mesh.compute_incident_entities(domain.topology, cell_outside, domain.topology.dim, domain.topology.dim-1)
left_facets= dolfinx.mesh.compute_incident_entities(domain.topology, cell_left, domain.topology.dim, domain.topology.dim-1)
common_facets = np.intersect1d(inside_facets, outside_facets)

integration_entities=internal_side_order(domain,common_facets,ct.values)

nfacets=domain.topology.index_map(fdim).size_local+domain.topology.index_map(fdim).num_ghosts
remaining_facets=np.arange(nfacets, dtype=np.int32)
remaining_facets_logic=np.full_like(remaining_facets,True,dtype=np.bool_)
remaining_facets_logic[dolfinx.mesh.exterior_facet_indices(domain.topology)]=False
remaining_facets_logic[common_facets]=False
remaining_facets=remaining_facets[remaining_facets_logic]
integration_entities2=internal_side_order(domain,remaining_facets,ct.values)

dS = ufl.Measure("dS", domain=domain,  subdomain_data=[[1,np.array(integration_entities2, dtype=np.int32)],[2,np.array(integration_entities, dtype=np.int32)]])
dx = ufl.dx

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = dolfinx.fem.Constant(domain, dolfinx.default_scalar_type(0))
leftValue=dolfinx.fem.Constant(domain, dolfinx.default_scalar_type(-1))
rightValue=dolfinx.fem.Constant(domain, dolfinx.default_scalar_type(1))

n=ufl.FacetNormal(domain)
h = ufl.Circumradius(domain)
alpha = dolfinx.fem.Constant(domain, dolfinx.default_scalar_type(100))
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * dx   

a += -ufl.inner(ufl.jump(v,n),ufl.avg(ufl.grad(u)))*dS(1)
a += -ufl.inner(ufl.avg(ufl.grad(v)),ufl.jump(u,n))*dS(1)
a +=  alpha/ufl.avg(h)*ufl.inner(ufl.jump(v,n),ufl.jump(u,n)) * dS(1)

a += - ufl.inner(n('+'), ufl.grad(u('+'))) * v('+') * (dS(2))
a += - ufl.inner(n('+'), ufl.grad(v('+'))) * u('+') * (dS(2)) + alpha/h('+')  * ufl.inner(u('+'), v('+')) * (dS(2))

a += - ufl.inner(n('-'), ufl.grad(u('-'))) * v('-') * (dS(2))
a += - ufl.inner(n('-'), ufl.grad(v('-'))) * u('-') * (dS(2)) + alpha/h('-')  * ufl.inner(u('-'), v('-')) * (dS(2))

L = f * v * dx
L += - ufl.inner(n('+'), ufl.grad(v('+'))) * leftValue * dS(2) + alpha/h('+')  * ufl.inner(leftValue, v('+')) * dS(2)
L += - ufl.inner(n('-'), ufl.grad(v('-'))) * rightValue * dS(2) + alpha/h('-')  * ufl.inner(rightValue, v('-')) * dS(2)

A = dolfinx.fem.petsc.create_matrix(dolfinx.fem.form(a))
b = dolfinx.fem.petsc.create_vector(dolfinx.fem.form(L))
dolfinx.fem.petsc.assemble_vector(b, dolfinx.fem.form(L))
dolfinx.fem.petsc.assemble_matrix(A, dolfinx.fem.form(a))
A.assemble()
ksp=PETSc.KSP()
solver = ksp.create(domain.comm)
solver.setOperators(A)
uhd = dolfinx.fem.Function(V)
solver.solve(b, uhd.x.petsc_vec)
5

Remaining domains meaning all other internal surfaces (those not marked with 8 and 10)?

Instead of creating a constant on each side of the surface, I would create a DG-0 function, and interpolate into the cells on that side the constant you would like to apply:

first_cells = integration_entities.reshape(-1, 4)[0]
second_cells = integration_entities.reshape(-1, 4)[2]
Q = dolfinx.fem.functionspace(mesh, ("DG", 0))
q.interpolate(lambda x: np.full(x.shape[1], -1), cells=first_cells))
q.interpolate(lambda x: np.full(x.shape[1], 1), cells=second_cells))
q.x.scatter_forward()

Then you should be able to write

as
L+= ufl.jump(ufl.inner(n, ufl.grad(v))*q+ alpha/h*q*v)*dS(2)

6

I mean the interior facets that in your case is not marked 12.

I have done it in the code above by taking all facets, deleting the exterior facets and the common_facets. Maybe there is a simpler way?

I understand now your dg-0 comment, Thanks

/Torben