How to enumerate the cells of a mesh with dolfinx v0.6.0?

General
1

I have an example I am cobbling together:

import numpy as np
import ufl
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from ufl import ds, dx, grad, inner
from dolfinx import fem, io, mesh, plot

msh = mesh.create_box(MPI.COMM_WORLD, [np.array([0.0, 0.0, 0.0]),
                                  np.array([2.0, 1.0, 1.0])], [16, 16, 16],
                 mesh.CellType.tetrahedron)

V = fem.FunctionSpace(msh, ("Lagrange", 1))


facets = mesh.locate_entities_boundary(msh, dim=(msh.topology.dim - 1),
                                       marker=lambda x: np.logical_or(np.isclose(x[0], 0.0),
                                                                      np.isclose(x[0], 2.0)))

dofs = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facets)

bc = fem.dirichletbc(value=ScalarType(0), dofs=dofs, V=V)

#for i, cell in enumerate(cells(mesh)):
#    if i==1: # for the second element
#        A = assemble_local(inner(grad(u), grad(v)) *dx, cell)

I got the comment prefixed code to enumerate cells from another post of the fenics discourse forum. It seems that cells and assemble_local are no longer accessible by those names.

I scanned the docker source for assemble_local and didn’t find anything there any longer. The cells function seems to be found in mesh, geometry and fem modules.

Is there a replacement in latest docker source for assemble_local someplace (maybe like some version of assemble…)?

Is fem.cells a suitable way to replace the comment prefixed code I found it to enumerate each cell?

(Wasn’t really officially part of my question but I feel curious here…),
(In terms of the usage of locate_entities_boundary that came from poisson example it seems that the lambda is only looking to the X co-ordinate for a bounds checking. I feel perplexed how that works without accounting for Y. Any insight on this? Suppose a box that lambda is moving to from a rectangle should checking just X with the lambda be sufficient for use with dolfinx or will other dimensions need to be looked at in terms of this style of coding with box now?)

2

It is very simply to enumerate cells in DOLFINx. The number of cells on a process can be found by
num_cells = msh.topology.index_map(msh.topology.dim).size_local, where the ith local index is i.

It is not included in DOLFINx by default. However, here is a code that does what you want

import numpy as np
import ufl
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from ufl import ds, dx, grad, inner
from dolfinx import fem, io, mesh, plot
import cffi



class LocalAssembler():
    def __init__(self, form):
        self.form = fem.form(form)
        self.update_coefficients()
        self.update_constants()
        subdomain_ids = self.form.integral_ids(fem.IntegralType.cell)
        assert(len(subdomain_ids) == 1)
        assert(subdomain_ids[0] == -1)
        is_complex = np.issubdtype(ScalarType, np.complexfloating)
        nptype = "complex128" if is_complex else "float64"
        self.kernel = getattr(self.form.ufcx_form.integrals(fem.IntegralType.cell)[0], f"tabulate_tensor_{nptype}")
        self.active_cells = self.form.domains(fem.IntegralType.cell, -1)
        assert len(self.form.function_spaces) == 2
        self.local_shape = [0,0]
        for i, V in enumerate(self.form.function_spaces):
            self.local_shape[i] = V.dofmap.dof_layout.block_size * V.dofmap.dof_layout.num_dofs

        e0 = self.form.function_spaces[0].element
        e1 = self.form.function_spaces[1].element
        needs_transformation_data = e0.needs_dof_transformations or e1.needs_dof_transformations or \
            self.form.needs_facet_permutations
        if needs_transformation_data:
            raise NotImplementedError("Dof transformations not implemented")

        self.ffi = cffi.FFI()
        V = self.form.function_spaces[0]
        self.x_dofs = V.mesh.geometry.dofmap

    def update_coefficients(self):
        self.coeffs = fem.assemble.pack_coefficients(self.form)[(fem.IntegralType.cell, -1)]

    def update_constants(self):
        self.consts = fem.assemble.pack_constants(self.form)

    def update(self):
        self.update_coefficients()
        self.update_constants()

    def assemble_matrix(self, i:int):

        x = self.form.function_spaces[0].mesh.geometry.x
        x_dofs = self.x_dofs.links(i)
        geometry = np.zeros((len(x_dofs), 3), dtype=np.float64)
        geometry[:, :] = x[x_dofs]

        A_local = np.zeros((self.local_shape[0], self.local_shape[0]), dtype=ScalarType)
        facet_index = np.zeros(0, dtype=np.intc)
        facet_perm = np.zeros(0, dtype=np.uint8)
        if self.coeffs.shape == (0, 0):
            coeffs = np.zeros(0, dtype=ScalarType)
        else:
            coeffs = self.coeffs[i,:]
        ffi_fb = self.ffi.from_buffer
        self.kernel(ffi_fb(A_local), ffi_fb(coeffs), ffi_fb(self.consts), ffi_fb(geometry),
               ffi_fb(facet_index), ffi_fb(facet_perm))
        return A_local


msh = mesh.create_box(MPI.COMM_WORLD, [np.array([0.0, 0.0, 0.0]),
                                  np.array([2.0, 1.0, 1.0])], [16, 16, 16],
                 mesh.CellType.tetrahedron)

V = fem.FunctionSpace(msh, ("Lagrange", 1))

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

assembler = LocalAssembler(inner(grad(u), grad(v)) *dx)

for cell in range(msh.topology.index_map(msh.topology.dim).size_local):
    if cell == 2:
        A = assembler.assemble_matrix(2)
        print(A)

Might I ask what you want to use the local assembler for?

locate_entities_boundary takes in a two dimensional array of shape (3, num_points), thus x[0] is the x coordinate, x[1] the y coordinate and x[2] the z-coordinate.

4 Likes
3

So about the local assembler, so far I have Mark S. Gockenbach’s book (Understanding the Finite Element Method). So then, he talks about a Gaussian Quadrature for Triangle a 2D shape. So in my efforts to understand the 3D realm so far I found out that there also is one a Gaussian Quadrature rule set for tetrahedron so I am hoping that the local assembler will help me to understand how those things are working. (Mark S. also talks about an isoparametric approach)

So thank you sou much for asking because the next question I would have is what methodology is dolfinx using for tet in 3D (be it isoparametric, Gaussian Quadrutre rule set, Some other approach etc.). (Granted this in the source code to look at you know its fairly large set of code…)

(I do have Josep Birnic’s novel which does forumulate a K stifness matrix for tet in a much simpler to understand approach using L1,L2,L3,L4 values however at this time I don’t feel sure how Birnic’s approach can tie into dolfinx usage.)

4

FEniCSx uses iso-parametric elements, meaning that you can describe the cell with a P-th order Lagrange element (usually with equispaced dofs).

The quadrature scheme is chosen depending on the form (you can also set these yourself, see: Integration of matrix elements - #2 by dokken and its referenced) and the integrals are performed by mapping the integral to the reference element.

If you want to go into how the code actually looks like, you can formulate everything with ufl and Basix, as done in: dolfinx/cpp/demo/poisson/poisson.py at main · FEniCS/dolfinx · GitHub
and call python3 -m ffcx poisson.py to generate a header file with C code for assembling the local tensor of a and L

1 Like
5

Hello, Mr. dokken. I want to get the element stiffness matrix of the specify cell using dolfinx0.7.2. And I run the code you had post above. But it out put the following error:

Traceback (most recent call last):
  File "/mnt/d/BaiduNetdiskWorkspace/MultiscaleDynamicTOP/code/FEniCSBasedProject/Load_code/ke.py", line 78, in <module>
    assembler = LocalAssembler(inner(grad(u), grad(v)) *dx)
  File "/mnt/d/BaiduNetdiskWorkspace/MultiscaleDynamicTOP/code/FEniCSBasedProject/Load_code/ke.py", line 14, in __init__
    self.update_coefficients()
  File "/mnt/d/BaiduNetdiskWorkspace/MultiscaleDynamicTOP/code/FEniCSBasedProject/Load_code/ke.py", line 40, in update_coefficients
    self.coeffs = fem.assemble.pack_coefficients(self.form)[(fem.IntegralType.cell, -1)]
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/assemble.py", line 81, in pack_coefficients
    return _pack(form)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/assemble.py", line 79, in _pack
    return _pack_coefficients(form)
TypeError: pack_coefficients(): incompatible function arguments. The following argument types are supported:
    1. (form: dolfinx::fem::Form<float, float>) -> Dict[Tuple[dolfinx::fem::IntegralType, int], numpy.ndarray[numpy.float32]]
    2. (form: dolfinx::fem::Form<double, double>) -> Dict[Tuple[dolfinx::fem::IntegralType, int], numpy.ndarray[numpy.float64]]
    3. (form: dolfinx::fem::Form<std::complex<float>, float>) -> Dict[Tuple[dolfinx::fem::IntegralType, int], numpy.ndarray[numpy.complex64]]
    4. (form: dolfinx::fem::Form<std::complex<double>, double>) -> Dict[Tuple[dolfinx::fem::IntegralType, int], numpy.ndarray[numpy.complex128]]

Invoked with: <dolfinx.fem.forms.Form object at 0x7f2c21bbb3a0>

Did you forget to `#include <pybind11/stl.h>`? Or <pybind11/complex.h>,
<pybind11/functional.h>, <pybind11/chrono.h>, etc. Some automatic
conversions are optional and require extra headers to be included
when compiling your pybind11 module.

Process finished with exit code 1

Can you tell me what’s going on here?

6

@dokken Excuse me, Mr. dokken. Can you help me see what is wrong with this. I’m trying to get the cell stiffness matrix using dolfinx 0.7.3, but it’s reporting an error.

7

Please do not tag me in posts. I try to answer as many posts on this forum as I have time for.
Here is a slightly updated version (for 0.7.0),
but it includes assembly over exterior facets as an option

# Local assembly of cell integral or exterior facet integrals
# Author: Jørgen S. Dokken
# SPDX-License-Identifier: MIT

import numpy as np
import ufl
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from ufl import  dx, grad, inner, ds
from dolfinx import fem, mesh
import cffi


_type_to_offset_index = {fem.IntegralType.cell: 0, fem.IntegralType.exterior_facet: 1}


class LocalAssembler():

    def __init__(self, form, integral_type: fem.IntegralType =fem.IntegralType.cell):
        self.form = fem.form(form)
        self.integral_type = integral_type
        self.update_coefficients()
        self.update_constants()

        subdomain_ids = self.form._cpp_object.integral_ids(integral_type)
        assert(len(subdomain_ids) == 1)
        assert(subdomain_ids[0] == -1)
        is_complex = np.issubdtype(ScalarType, np.complexfloating)
        nptype = "complex128" if is_complex else "float64"
        o_s = self.form.ufcx_form.form_integral_offsets[_type_to_offset_index[integral_type]]
        o_e = self.form.ufcx_form.form_integral_offsets[_type_to_offset_index[integral_type]+1]
        assert o_e - o_s == 1

        self.kernel = getattr(self.form.ufcx_form.form_integrals[o_s], f"tabulate_tensor_{nptype}")
        self.active_cells = self.form._cpp_object.domains(integral_type, -1)
        assert len(self.form.function_spaces) == 2
        self.local_shape = [0,0]
        for i, V in enumerate(self.form.function_spaces):
            self.local_shape[i] = V.dofmap.dof_layout.block_size * V.dofmap.dof_layout.num_dofs

        e0 = self.form.function_spaces[0].element
        e1 = self.form.function_spaces[1].element
        needs_transformation_data = e0.needs_dof_transformations or e1.needs_dof_transformations or \
            self.form._cpp_object.needs_facet_permutations
        if needs_transformation_data:
            raise NotImplementedError("Dof transformations not implemented")

        self.ffi = cffi.FFI()
        V = self.form.function_spaces[0]
        self.x_dofs = V.mesh.geometry.dofmap

    def update_coefficients(self):
        self.coeffs = fem.assemble.pack_coefficients(self.form._cpp_object)[(self.integral_type, -1)]

    def update_constants(self):
        self.consts = fem.assemble.pack_constants(self.form._cpp_object)

    def update(self):
        self.update_coefficients()
        self.update_constants()

    def assemble_matrix(self, i:int, local_index:int = 0):

        x = self.form.function_spaces[0].mesh.geometry.x
        x_dofs = self.x_dofs[i]
        geometry = np.zeros((len(x_dofs), 3), dtype=np.float64)
        geometry[:, :] = x[x_dofs]

        A_local = np.zeros((self.local_shape[0], self.local_shape[0]), dtype=ScalarType)
        facet_index = np.array([local_index], dtype=np.intc)
        facet_perm = np.zeros(0, dtype=np.uint8)
        if self.coeffs.shape == (0, 0):
            coeffs = np.zeros(0, dtype=ScalarType)
        else:
            coeffs = self.coeffs[i,:]
        ffi_fb = self.ffi.from_buffer
        self.kernel(ffi_fb(A_local), ffi_fb(coeffs), ffi_fb(self.consts), ffi_fb(geometry),
               ffi_fb(facet_index), ffi_fb(facet_perm))
        return A_local


msh = mesh.create_unit_square(MPI.COMM_WORLD, 3, 3,
                 mesh.CellType.triangle)

V = fem.FunctionSpace(msh, ("Lagrange", 1))

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

assembler = LocalAssembler(inner(grad(u), grad(v)) *dx)

for cell in range(msh.topology.index_map(msh.topology.dim).size_local):
    if cell == 2:
        A = assembler.assemble_matrix(2)
        print(A)



facet_assembler = LocalAssembler(inner(u, v) *ds, fem.IntegralType.exterior_facet)

for cell in range(msh.topology.index_map(msh.topology.dim).size_local):
    if cell == 2:
        A = facet_assembler.assemble_matrix(2, local_index = 0)
        print(A)

yielding

[[ 0.5 -0.5  0. ]
 [-0.5  1.  -0.5]
 [ 0.  -0.5  0.5]]
[[0.         0.         0.        ]
 [0.         0.11111111 0.05555556]
 [0.         0.05555556 0.11111111]]
8

Thank you very much. And I am sorry for what I have done.