notes on FEniCSx

Installation

  • Currently developed actively, and api also changes
  • has complex mode and real mode; PDE solver library PETSc has either real or complex build, and FEniCSx has to be switched to the corresponding mode. Below command to check the available version and build, and install PETSc(Corresponding FEniCSx will be automatically installed). It is an idea to install both in different environment to allow switching.

’’’ conda install -c conda-forge fenics-dolfinx==0.9.0 mpich pyvista conda search -c conda-forge PETSc conda install -c conda-forge petsc=3.22.3=complex_h1434fe3_0 (for 0.8.0, conda install -c conda-forge petsc=3.21.6=complex_h9a830cf_0) conda install -c conda-forge meshio conda install -c conda-forge jupyterlab h5py jupyterlab_widgets conda install -c conda-forge matplotlib ipywidgets plotly cd viskex python3 -m pip install ‘.[tutorials]’ source /opt/homebrew/Caskroom/miniconda/base/envs/FEMBEM25Mar which python pip install bempp-cl ‘’’

Coding note

tabulation

evaluate the basis function at a certain set of points in the reference element by using tabulation. The first argument: the number of spatial derivatives of the basis functions we want to compute The second argument: a set of input points to evaluate the basis functions at as a numpy array of shape (num_points, reference_cell_dimension)

points = np.array([[0.0, 0.1], [0.3, 0.2]])
values = element.tabulate(0, points)

Outputs the array of shape (num_spatial_derivatives+1, num_points, num_basis_functions)

In some cases, one might want to use a lower accuracy for tabulation of basis functions to speed up computations. This can be changed in basix by adding dtype=np.float32 to the element constructor. Basix also allows other extra options to tweak finite elements–try changing the options and see the shape of the shape functions

Mapping

As we have already seen, we can describe any cell in our subdivided domain with a mapping from the reference element. However, as we want to integrate over each element individually, we need to map the basis functions to and from the reference element. We call this map: .

What property does the covariant Piola mapping preserve?
We observe that the covariant Piola map maps normals to normals, since a 0 tangential component is mapped to a 0 tangent.

Entity map

This guides you how to create entity map to assemble across subdomain and main domain. Entity map is the reverse dictionary(sub->parent domain) of the parent->sub mapping obtained as the 2nd return of mesh.create_submesh.

creating entities and connectivity/ obtaining boundary facets

mesh.topology.create_entities(mesh.topology.dim - 1)
gmshmesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)

When a mesh is created in DOLFINx we only compute the unique global numbering and ownership of cells and vertices. However, one can compute this numbering for any sub-entity by calling mesh.topology.create_entities(j). For instance computing the ownership and numbering of the facets can be done with mesh.topology.create_entities(mesh.topology.dim-1).

When a mesh is created in DOLFINx, the relation between the cells and the vertices is computed. However, one can compute the relationship between any set of entities in the mesh with dolfinx.mesh.create_connectivity(i, j) where i is the dimension of the index to map from and j is the dimension of the entities to map from. As an example, calling mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim) will compute the relationship between all facets in the mesh and their cells. The inverse map (cell-to-facet) is computed with mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim)-1

With these two preliminaries computed, DOLFINx can determine which facets that are connected to a single cell only.

## Another feature of variational formulations is that the test function v is required to vanish on the parts of the boundary where the solution u is known. See for instance [LM19].

This simple but very powerful method for constructing test problems is called the method of manufactured solutions. First pick a simple expression for the exact solution, plug into the equation to obtain the right-hand side (source term f ). Then solve the equation with this right hand side, and using the exact solution as boundary condition. Finally, we create a program that tries to reproduce the exact solution.

Entities vs degrees-of-freedom Entities: A mesh entity in FEniCS is either a vertex, an edge, a face, or a cell (triangle or tetrahedron) Degrees-of-freedom: mesh nodes and interpolation points ct vs ft ct: volumetric element eg)hexahedron ft: surface element eg)quadrilateral

Creating a source on a pointO

Next, we create a function containing the exact solution (which will also be used in the Dirichlet boundary condition) and the corresponding source function for the right hand side. Note that we use ufl.SpatialCoordinate to define the exact solution, which in turn is interpolated into uD and used to create the source function

Defining a boundary condition for mixed elements

V.sub(0).sub(0) allows accessing the item of a list that’s a part of the mixed element.

    clamped_dofsW = fem.locate_dofs_topological(V.sub(0), facet_dim, clamped_facets)
    bcs = [fem.dirichletbc(PETSc.ScalarType(0.0), clamped_dofsW, V.sub(0))]
    clamped_dofsT = fem.locate_dofs_topological(V.sub(1), facet_dim, clamped_facets)
    bcs += [fem.dirichletbc(PETSc.ScalarType(0.0), clamped_dofsT, V.sub(1).sub(0))]
    bcs += [fem.dirichletbc(PETSc.ScalarType(0.0), clamped_dofsT, V.sub(1).sub(1))]

Avoid mixed elements for coupled PDEs

When assembling the matrix blocks, mixed elements block is created to construct its own diagonals, ending up with a different order of dofs–BCs are applied to wrong dofs in the end.

When porting codes (from FEniCSx version, from different python environment)

When modifying the code written for real mode to run on complex mode

- Wrap constant with PETSc.ScalarType [Ref](https://jsdokken.com/dolfinx-tutorial/chapter1/complex_mode.html#variational-problem/)

	rho = fem.Constant(domain, PETSc.ScalarType(2700.0))

- Expected error [ incompatible function arguments. The following argument types are supported: ]
- conjugate tensor product either by ufl.conj or ufl.inner(second operand conjugated)	
	- Expected error [ ArityMismatch: Failure to conjugate test function in complex Form ]
- [# Similarly, if we want to use the function `ufl.derivative` to take derivatives of functionals, we need to take some special care. As `ufl.derivative` inserts a `ufl.TestFunction` to represent the variation, we need to take the conjugate of this to be able to use it to assemble vectors.](https://github.com/jorgensd/dolfinx-tutorial/blob/main/chapter1/complex_mode.py)
	
	
- Defining boundary condition with constants, below failed in solving the problem due to an incompatible type error

bcs = [
    fem.dirichletbc(np.zeros((2,)), left_dofs, V),
    fem.dirichletbc(np.zeros((2,)), right_dofs, V),
]

- Changing to numpy to PETSc.ScalarType resolved the issue

bcs = [
    fem.dirichletbc(PETSc.ScalarType([0, 0]), left_dofs, V),
    fem.dirichletbc(PETSc.ScalarType([0, 0]), right_dofs, V),
]

PETSc.ScalarType wraps the constant source on the right hand side. This is because we want the integration kernels to assemble into the correct floating type.

  • When plotting, some values have to be set back to real values before plot
grid.point_data["Deflection"] = u.x.array.reshape(-1, 3).real

When modifying from one version to another

  • ERROR [compute_integration_domains(): incompatible function arguments.]
    • Some versions needs explicit type definition for meshtags
      facet_tags = meshtags(mesh, mesh.topology.dim - 1,
np.concatenate([facets_left]),
    np.array([0] * len(facets_left)).astype('int32'))

FEniCSx API changne:Misc

FEniCSx API change:Function

When getting the value of mesh

display(kappa.x.array[:])

When importing mesh to FEniCSx

  • gmshio seems to fail to read hexahedron…Unknown cell type-KeyError: np.int32(6)
    • works fine with tetrahedron
  • cell:Volume physical group number

  • facet_tags: Surface physical number

  • viskex can only plot dolfinx.mesh.Mesh object, not meshio mesh object
    • so better to write in xdmf and import from dolfinx <- not necessaryly cf)dolfinx.io.gmshio.model_to_mesh
    • also, this process helps separating different elemtns of the mesh, such as hexa, triangle, quad and line
    • Importantly, xdmf can not handle mixed elements, so different elements have to be stored separately

Setting boundary condition

def Omega_0(x):
    return x[1] <= 0.5


def Omega_1(x):
    return x[1] >= 0.5

Note that both functions use a ≤ or ≥, as FEniCSx will evaluate each cell at all of the vertices, and thus has to return True for all vertices aligned with the interface to be marked properly.

Defining material property

  • when density is a part of the sound source term, the entire volume domain will be integrated without density in the integrand–move the denisty to the LHS so that density can be defined as a part of volume integration where different density can be defined on different part of the volume.

Dealing with Mesh

Check the shape of the imported mesh

msh, cell, facet_tags=gmshio.read_from_msh('./data/MembraneAbsorber.msh',MPI.COMM_WORLD,0,gdim=3)
topology, cell_types, geometry = plot.vtk_mesh(msh, 3)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
p = pyvista.Plotter()
p.add_mesh(grid, show_edges=True)
p.view_xy()
p.show_axes()
p.show()

Check the size of the imported mesh

create measures for surface and volume, and then integrate over the region to calculate the area and volume–making sure the size of the mesu is correct

msh, cell, facet_tags=gmshio.read_from_msh('./data/MembraneAbsorber.msh',MPI.COMM_WORLD,0,gdim=3)
display(cell.values,facet_tags.values)
ds=Measure('ds',domain=msh,subdomain_data=facet_tags)
dx=Measure('dx',domain=msh,subdomain_data=cell)
AreaTop=fem.assemble_scalar(fem.form(1.0*ds(2)))
VPoro=fem.assemble_scalar(fem.form(1.0*dx(1)))
VAir=fem.assemble_scalar(fem.form(1.0*dx(100)))

display('Expecting:0.8x0.8=0.064:',AreaTop)
display('Expecting:0.05x0.8:0.04',VPoro)
display('Expecting:0.75x0.8:0.6',VAir)

Solving a coupled PDEs

Monolithic Direct Solver

  • Assemble all the unknowns in a single large matrix to be solved simulatneously using a direct solver such as LU decomposition
  • simple to implement, free from convergence issues
  • Expensive for large problems, memory intensive

Nested (block) matrix solver

  • Split the matrix into blocks and solve each separately and iteratively
  • One unknown is solved while other unknowns are treated as preconditioners?
  • Memory efficient, more suitable for large problems, specialized solvers can be employed for each block
  • Complex implementation, could be slower for small problems

Iterative solver (Schur complement)

  • Instead of solving the entire system directly, solves for one unknown at a time using Schur complement
  • Highly efficient and scalable for larg systems
  • Complex implementation and the performance is reliant on the choice of preconditioners

Non-block direct solver

  • Solves like monolithic block direct solver, but the structure of the system is not preserved and is treated as a single dense system
  • easy to implement, as the whole matrix is solved with a single solver
  • less efficient for large problem

Zakki

Jupyter

COMSOL

  • COMSOL6.3-Shell and Plate introduction-
    The Shell interface ( ) is intended for mechanical analysis of thin-walled structures. The formulation used in the Shell interface is a Mindlin–Reissner type, which means that transverse shear deformations are accounted for. It can be used for rather thick shells as well as thin ones. It is possible to prescribe an offset in a direction normal to a selected surface, which for example can be used when meshing imported geometries. The Shell interface also includes other features such as damping, thermal expansion, and initial stresses and strains.

  • Structural Mechanics Theory(p.494 for elastic moduli)

  • Equation View

FEniCSx

Discourse

Peripherals

-DefElement:an encyclopedia for finite element definitions

Bempp

Physics

Misc

V_{sphere} = \frac{4}{3}\pi r^3

Everything else

Studio design

Acoustics

Manufactures

Composers/Sound designers

Rhino

Utility softwares and libraries

Web template

Shelves