elastodynamicsx.utils

Module contents

The utils module contains various tools that do not fit into the other packages

elastodynamicsx.utils.make_facet_tags(domain: Mesh, boundaries: Tuple[Tuple[int, Callable]]) → MeshTags[source]: Shortcut for make_tags(domain, locators, type_=’boundaries’)

elastodynamicsx.utils.make_cell_tags(domain: Mesh, subdomains: Tuple[Tuple[int, Callable]]) → MeshTags[source]: Shortcut for make_tags(domain, locators, type_=’domains’)

elastodynamicsx.utils.make_tags(domain: Mesh, locators: Tuple[Tuple[int, Callable]], type_='unknown') → MeshTags[source]

Parameters:

domain – A mesh
locators – A tuple of tuples of the type (tag, fn), where: tag: is an int to be assigned to the cells or facets fn: is a boolean function that take ‘x’ as argument
type – either ‘domains’ or ‘boundaries’

Returns:

A MeshTags object (dolfinx.mesh)

Adapted from:: https://jsdokken.com/dolfinx-tutorial/chapter3/robin_neumann_dirichlet.html

Example

import numpy as np
from mpi4py import MPI
from dolfinx.mesh import create_unit_square
#
domain = create_unit_square(MPI.COMM_WORLD, 10, 10)

boundaries = ((1, lambda x: np.isclose(x[0], 0)),
              (2, lambda x: np.isclose(x[0], 1)),
              (3, lambda x: np.isclose(x[1], 0)),
              (4, lambda x: np.isclose(x[1], 1)))
facet_tags = make_tags(domain, boundaries, 'boundaries')

Omegas = ((1, lambda x: x[1] <= 0.5),
          (2, lambda x: x[1] >= 0.5))
cell_tags = make_tags(domain, Omegas, 'domains')

class elastodynamicsx.utils.ParallelEvaluator(domain: Mesh, points: ndarray, padding: float = 0.0001)[source]

Bases: object

Convenience class to evaluate functions (fem.Function) when the mesh is distributed over several processes

Parameters:

domain – a distributed mesh (MPI.COMM_WORLD)
points – shape=(3, nbpts). Should be defined either for all processes, or for the proc no 0. Does not support ‘points’ being scattered on several processes.
padding – close-to-zero parameter used in dolfinx.cpp.geometry.determine_point_ownership

Example

# Define some output points
x = np.linspace(0, 1, num=30)
y = np.zeros_like(x)
z = np.zeros_like(x)
points = np.array([x, y, z])

# Initialize the evaluator
paraEval = ParallelEvaluator(domain, points)

# Perform function evaluations
u_eval_local = u.eval(paraEval.points_local, paraEval.cells_local)

# Gather all procs
u_eval = paraEval.gather(u_eval_local, root=0)

# Do something, e.g. export to file
if domain.comm.rank == 0:
    np.savez('u_eval.npz', x=points.T, u=u_eval)

Adapted from:: https://github.com/jorgensd/dolfinx-tutorial/issues/116

property nb_points_local: int

gather(eval_results: ndarray, root=0) → ndarray | None[source]

elastodynamicsx.utils.spectral_element(name: str, cell_type, degree: int, shape: Tuple[int, ...] | None = None) → _BasixElement[source]

A spectral element that can be used in a dolfinx.fem.FunctionSpace

Parameters:

name – One of (“GLL”, “GL”, “Legendre”)
cell_type – Elements can be defined from any type, but diagonal mass matrices can only be obtained using GLL / GL quadratures that require cell types interval (1D) / quadrilateral (2D) / hexahedron (3D)
degree – The maximum degree of basis functions

Example

from mpi4py import MPI
from dolfinx import mesh, fem
from elastodynamicsx.utils import spectral_element, spectral_quadrature
#
degree = 4
specFE = spectral_element("GLL", mesh.CellType.quadrilateral, degree)
specmd = spectral_quadrature("GLL", degree)
domain = mesh.create_unit_square(MPI.COMM_WORLD, 10, 10,
                                 cell_type=mesh.CellType.quadrilateral)
V = fem.functionspace(domain, specFE)

#######
# Compile mass matrices using pure fenicsx code
import ufl
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)

# here we do not specify quadrature information -> non-diagonal mass matrix
a1_non_diag = fem.form(ufl.inner(u,v) * ufl.dx)
M1_non_diag = fem.petsc.assemble_matrix(a1_non_diag)
M1_non_diag.assemble()

# here we specify the quadrature metadata -> the matrix becomes diagonal
a2_diag     = fem.form(ufl.inner(u,v) * ufl.dx(metadata=specmd))
M2_diag     = fem.petsc.assemble_matrix(a2_diag)
M2_diag.assemble()

#######
# Compile mass matrices using elastodynamicsx.pde classes
from elastodynamicsx.pde import PDE, material

# here we do not specify quadrature information -> non-diagonal mass matrix
mat3          = material(V, 'scalar', 1, 1)
pde3_non_diag = PDE(V, [mat3])
M3_non_diag   = pde3_non_diag.M()

# elastodynamicsx offers two ways to specify the quadrature metadata

# 1. by passing the metadata as kwargs to each material / boundary condition / bodyforce

mat4          = material(V, 'scalar', 1, 1, metadata=specmd)
pde4_diag     = PDE(V, [mat4])
M4_diag       = pde4_diag.M()

# 2. by passing the metadata to the PDE.default_metadata, which is being used by
# all materials / boundary conditions / bodyforces at __init__ call

from elastodynamicsx.pde import PDECONFIG
PDECONFIG.default_metadata = specmd  # default for all forms in the pde package (materials, BCs, ...)

mat5          = material(V, 'scalar', 1, 1)
pde5_diag     = PDE(V, [mat5])
M5_diag       = pde5_diag.M()

# spy mass matrices
from elastodynamicsx.plot import spy_petscMatrix
import matplotlib.pyplot as plt
for i,M in enumerate([M1_non_diag, M2_diag, M3_non_diag, M4_diag, M5_diag]):
    fig = plt.figure()
    fig.suptitle('M'+str(i+1))
    spy_petscMatrix(M)
plt.show()

elastodynamicsx.utils.spectral_quadrature(name: str, degree: int) → dict[source]

A quadrature metadata to build diagonal mass matrices when used with corresponding spectral elements.

Parameters:

name – One of (“GLL”, “GL”, “Legendre”)
degree – The maximum degree of basis functions

Example

See doc of ‘spectral_element’.