Incorrect result for assembly over boundary part

Errors
1

Hi,

I’m having trouble with the following code, specifically when I want to take the integral over a part of the boundary. The solver part of the code works and I can view the results as expected in ParaView. I’d like to take the average of the field u over several parts of the boundary (indexed 3 - 18). However, I get the following result when I run the code snippet below, which doesn’t match up with the ParaView visualization. This seems like I’m probably misunderstanding something about the assemble function. I’m not sure it does, but if it helps, print(u*ds(5)) results in “{ f_25 } * ds(<Mesh #0>[5], {})”.

Thanks in advance!

# %% Define Materials on Subdomains
class materials(UserExpression):
    """
    Assign a value to a parameter piece-wise to each subdomain
    """
    def __init__(self, subdomains, param_value_list, **kwargs):
        super().__init__(**kwargs)
        self.subdomains = subdomains
        self.val_dict = self.generate_dictionary(param_value_list)

    def generate_dictionary(self, param_values):
        val_dict = {}
        for sd,val in param_values:
            val_dict[f'{sd}'] = val
        return val_dict

    def eval_cell(self, values, x, cell):
        subdomain_label = self.subdomains[cell.index]
        values[0] = self.val_dict[f'{subdomain_label}']

# %% Import Mesh
output_directory='bundle_mesh'
msh, boundary, subdomain = load_xdmf(f'{output_directory}/mesh.xdmf',
                                      f'{output_directory}/boundary.xdmf', 
                                      f'{output_directory}/subdomain.xdmf')

# %% Define element
CG_elem = FiniteElement('CG', msh.ufl_cell(), 1)
V = FunctionSpace(msh, CG_elem)

dx = Measure('dx', domain=msh, subdomain_data=subdomain) 
ds = Measure('ds', domain=msh, subdomain_data=boundary) 

k = interpolate(materials(subdomain, [[1, Constant(1.0)], [2, Constant(1e6)]], element=CG_elem), V)

# %% Define Dirichlet Boundary Conditions
bcs = [DirichletBC(V, Constant(0.0), boundary, 7)] # Possible non-zero BC locations are indexed 3 through 18
g = -1.0

# %% Define Electrostatics Problem
u = TrialFunction(V)
v = TestFunction(V)
F = inner(k*grad(v), grad(u))*dx
F += v*g*ds(15) 

a = lhs(F)
L = rhs(F)

u = Function(V)
solve(a == L, u, bcs)
File(f'Simulation Results/u.pvd') << u

for boundary_part in np.unique(boundary.array()):
    facet_area = assemble(Constant(1.0)*ds(boundary_part))
    integrated_voltage = assemble(u*ds(boundary_part))
    print(boundary_part, facet_area, integrated_voltage)

Output:

‘’’
3 0.0 0.0
4 0.0 0.0
5 0.0 0.0
6 0.0 0.0
7 0.007362368556873695 0.0
8 0.0 0.0
9 0.0 0.0
10 0.0 0.0
11 0.0 0.0
12 0.0 0.0
13 0.0 0.0
14 0.0 0.0
15 0.0 0.0
16 0.0 0.0
17 0.0 0.0
18 0.0 0.0
‘’’

2

As your code is not complete and reproducable (several key definitions and the mesh is missing), it is very hard for people to assist you.

Please follow the guidelines in: Read before posting: How do I get my question answered?
and make sure that the issues you are getting are reproducable, i.e. that anyone can run the posted code and obtain the same results as you.

3

Thanks @dokken for your quick reply! I came across your the tutorials you created (The FEniCSx tutorial — FEniCSx tutorial) and I’ve been referencing them a lot lately–they’re super helpful! Thank you for those, too!

Regarding your response to my original comment, I put together a MWE below, but I was hoping the issue was related to the assembly of the integrals over the boundary_parts in the last for loop. I still don’t understand why the following code doesn’t output the length of the boundary_part labeled “7”.

assemble(1.0*ds(7))

What is the correct way to do this?

Thanks for your help! My MWE code is below:

from dolfin import *
import numpy as np

msh = UnitSquareMesh(10, 10)
class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1.0) and on_boundary

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0) and on_boundary

boundary = MeshFunction("size_t", msh, 1)
boundary.set_all(0)
Top().mark(boundary, 7)
Bottom().mark(boundary, 15)

dx = Measure('dx', domain=msh)
ds = Measure('ds', domain=msh, subdomain_data=boundary)

normal = FacetNormal(msh)

# %% Define element type
CG_elem = FiniteElement('CG', msh.ufl_cell(), 1)
V = FunctionSpace(msh, CG_elem)

k = Constant(1.0)

# %% Define Dirichlet Boundary Conditions
bcs = DirichletBC(V, 1.0, boundary, 7)
g = 1.0

# %% Define Electrostatics Problem
u = TrialFunction(V)
v = TestFunction(V)
a = k*inner(grad(u), grad(v))*dx
L = g*v*ds(15)

u = Function(V)
solve(a == L, u, bcs)
file = File(f'Simulation Results/u.pvd') 
file << u

for boundary_part in np.unique(boundary.array()):
    facet_area = assemble(1*ds(boundary_part))
    integrated_voltage = assemble(u*ds(boundary_part))
    print(boundary_part, facet_area, integrated_voltage)

import matplotlib.pyplot as plt
plot(u)
plt.show()

Output:
0 0.0 3.000000000000028
7 0.0 0.0
15 0.0 0.0

Figure:

test
test1920×1440 33.9 KB

4

This is a bug in FEniCS, due to the fact that the data in the boundary_parts is uint64 while dolfin uses ints.
A workaround is to call int(boundary_parts), as in the following code:

from dolfin import UnitSquareMesh, SubDomain, MeshFunction, Measure, assemble, near, Constant
import numpy as np

mesh = UnitSquareMesh(10, 10)
class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1.0) and on_boundary

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0) and on_boundary

boundary = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
boundary.set_all(0)
Top().mark(boundary, 1)
Bottom().mark(boundary, 2)

dx_c = Measure("dx", domain=mesh)
ds_c = Measure("ds", domain=mesh, subdomain_data=boundary)




for boundary_part in np.unique(boundary.array()):
    b_p = int(boundary_part)
    area = 1 * ds_c(b_p)
    print(b_p, assemble(area))
5

Thanks @dokken for your help on this. I tried your code and it does what I’m looking for, but if I change the domain labels from 1 and 2, such as below:

Top().mark(boundary, 7)
Bottom().mark(boundary, 15)

I end up with the output:
0 0.0
7 0.0
15 0.0

Is there something about the labels 1 and 2 for which the assembly works but not for other labels, such as 7 and 15? Do the boundary labels have to appear in sequential order (0, 1, 2, 3, …)? If this is the case, I’m a bit confused on what is happening on the back end because the integrals in the problem formulation seem be set up correctly because I can set Neumann boundary conditions as expected. If I build my mesh in pygmsh, do I have to go through the additional step of re-labeling the boundaries in a particular order? Can you provide guidance on this?

Thanks again!

6

You should remove the JIT cache, for instance by calling dijitso clean, as the integrals with 7 and 15 has already been compiled.

The order and number doesn’t matter as long as the cache is clean and you then call the forms with the int encapsulation around the marker value.

1 Like
7

Thank you for your help! I added the dijitso clean statement to my code as you suggested and it seems to have got it to work:

import os
os.system('dijitso clean')
for boundary_part in np.unique(boundary.array()):
    b_p = int(boundary_part)
    facet_area = assemble(1*ds(b_p))
    integrated_voltage = assemble(u*ds(b_p))
    print(boundary_part, facet_area, integrated_voltage)

That being said, is there a way to avoid doing this step if I set up my mesh/code better, or is this an inevitability if I want to pass information from these integrals over boundary parts to another bit of code?

Thanks again!

8

You should be able to remove thedijitso clean call now, as the cache has the correct integer input now.

I’m not sure if there is a pretty workaround (except doing int(marker).

EDIT : The code initally posted here did not work, as numpy casted the value to np.int64 as opposed to the native int. So the workaround is using int(marker).

2 Likes