Integration over a single cell or on a cell boundary

General
1

I would like to be able to integrate some quantities on a single cell and on cells boundaries, with the aim to compute residual error estimates. For that particular case, I know that the residuals can be computed by assembling a bilinear form on all the domain, however I want to try to compute them by hand, element by element.
So my question is how to compute integrals on a single cell or on two or three facets. By searching on Fenics literature and forums I learnt something, and I wrote the following code (which should compute the area of each cell and the internal boundary length of the cells, but, the code does not work):

from dolfin import *
from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt

squareg = Rectangle(Point(-1.0, -1.0), Point(1.0,1.0)) 
mesh2d = generate_mesh(squareg, 1)
plot(mesh2d)
plt.show()

area_cells = [0] * mesh2d.num_cells()
int_bound_length_cells = [0] * mesh2d.num_cells()
domain_markers = MeshFunction("bool", mesh2d,mesh2d.topology().dim())
boundary_markers = MeshFunction("bool",mesh2d,mesh2d.topology().dim()-1)
for cell in cells(mesh2d):
  domain_markers.set_all(False)
  boundary_markers.set_all(False)
  domain_markers[cell] = True
  for fct in facets(cell):
    if ~fct.exterior():  #### not working: facet connectivity not initialized
      boundary_markers[fct]=True  
  
  ## Just for testing ############
  c=0
  for cellj in cells(mesh2d):
    if domain_markers[cellj]: c+=1
  d=0
  for f in facets(mesh2d):
    if boundary_markers[f]: d+=1
  print("cell ", cell.index())
  print("cells integ: ", c)
  print("edg integ: ", d)
  ################################ 

  dx_cell = Measure('dx',domain=mesh2d, subdomain_data=domain_markers)
  ds_cell = Measure('ds',domain=mesh2d, subdomain_data=boundary_markers)
  area_cells[cell.index()] = assemble(1.0*dx_cell)  ##### gives error
  int_bound_length_cells[cell.index()] = assemble(1.0*ds_cell) #### gives error

print(area_cells)
print(int_bound_length_cells)

Running this code gives the following error message (which is due to the “assemble” command):

Traceback (most recent call last):
  File "/home/ariel/Documents/fenics/question_local_integral.py", line 38, in <module>
    area_cells[cell.index()] = assemble(1.0*dx_cell)  ##### gives error
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py", line 202, in assemble
    dolfin_form = _create_dolfin_form(form, form_compiler_parameters)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py", line 60, in _create_dolfin_form
    return Form(form,
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/form.py", line 82, in __init__
    self.set_cell_domains(subdomains)
TypeError: set_cell_domains(): incompatible function arguments. The following argument types are supported:
    1. (self: dolfin.cpp.fem.Form, arg0: dolfin.cpp.mesh.MeshFunctionSizet) -> None

Invoked with: <dolfin.fem.form.Form object at 0x7f3816287a60>, <dolfin.cpp.mesh.MeshFunctionBool object at 0x7f381e53e4b0>

I have two questions:

  1. I can not find the way to perform the local integrations. Is that possible? and how?
  2. (A minor problem) My code to check if a facet is on the exterior boundary of the domain (“if ~fct.exterior(): …”) is not working, and it seems that is because the “facet connectivity” is not initialized. How can that be done?

Thanks in advance.

2

Consider the following code (using a unit square, as I do not use mshr)

from dolfin import *
from fenics import *

mesh2d = UnitSquareMesh(10, 10)
area_cells = [0] * mesh2d.num_cells()
int_bound_length_cells = [0] * mesh2d.num_cells()
mesh2d.init(mesh2d.topology().dim()-1, mesh2d.topology().dim())
domain_markers = MeshFunction("size_t", mesh2d, mesh2d.topology().dim())
boundary_markers = MeshFunction("size_t", mesh2d, mesh2d.topology().dim()-1)
marker = 13
dx_cell = Measure('dx', domain=mesh2d,
                  subdomain_data=domain_markers, subdomain_id=marker)
ds_cell = Measure('ds', domain=mesh2d,
                  subdomain_data=boundary_markers, subdomain_id=marker)

for cell in cells(mesh2d):
    domain_markers.set_all(0)
    boundary_markers.set_all(0)
    domain_markers[cell] = marker
    for fct in facets(cell):
        if ~fct.exterior():
            boundary_markers[fct] = marker

    ## Just for testing ############
    c = 0
    for cellj in cells(mesh2d):
        if domain_markers[cellj]:
            c += 1
    d = 0
    for f in facets(mesh2d):
        if boundary_markers[f]:
            d += 1
    print("cell ", cell.index())
    print("cells integ: ", c)
    print("edg integ: ", d)
    ################################

    area_cells[cell.index()] = assemble(1.0*dx_cell)
    int_bound_length_cells[cell.index()] = assemble(1.0*ds_cell)

print(area_cells)
print(int_bound_length_cells)

It should resolve all your questions.

3

Thanks very much for your answer! I think I now understood how to define the Measure to use in assemble, and how to initialize the facet connectivity.
Sorry, that I have a again a question/remark. Running your code with a coarse mesh

mesh2d = UnitSquareMesh(2, 2)

with eight elements of area = 0.125, the first shown result is ok :

[0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125]

However the result about the lengths of the interior boundary of the cells gives:

[0.5, 0.5, 1.0, 0.0, 0.0, 1.0, 0.5, 0.5]

which is wrong since those boundaries have lengths (approx) 0.7071, 1.2071 or 1.7071.
It is strange also that the three edges are marked on each element (which is wrong, but furthermore does not corresponds with the result).
Reading your code I do not find inconsistencies, but the results are not what I expect.
Thanks again!

4

I missed the part saying you want to do interior integrals. If you want to integrate over interior facets, you should use the dS measure, not the ds measure.
See for instance Form language — Unified Form Language (UFL) 2021.1.0 documentation

5

Thank you!! That’s right