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

‘’’