elastodynamicsx.solvers

Module contents

The solvers module contains tools for solving PDEs

Time domain:

from elastodynamicsx.solvers import TimeStepper

Normal modes:

from elastodynamicsx.solvers import EigenmodesSolver

Frequency domain:

from elastodynamicsx.solvers import FrequencyDomainSolver

… in the future… Guided waves:

from elastodynamicsx.solvers import ...

class elastodynamicsx.solvers.TimeStepper(comm: Comm, tscheme: TimeScheme, **kwargs)[source]

Bases: object

Base class for solving time-dependent problems.

petsc_options_explicit_scheme = {'ksp_type': 'preonly', 'pc_type': 'lu'}

petsc_options_implicit_scheme_linear = {'ksp_type': 'preonly', 'pc_type': 'lu'}

build(**kwargs)[source]

Convenience static method that instanciates the required time-stepping scheme

Parameters:

args – (passed to the required scheme)

Keyword Arguments:

scheme – (required) Available options are ‘leapfrog’, ‘midpoint’ ‘linear-acceleration-method’, ‘newmark’, ‘hht-alpha’, ‘generalized-alpha’
**kwargs – (passed to the required scheme)

static Courant_number(domain: Mesh, c_max, dt)[source]

The Courant number: \(C = c_{max} \, \mathrm{d}t / h\), with \(h\) the cell diameter

Related to the Courant-Friedrichs-Lewy (CFL) condition.

Parameters:

domain – the mesh
c_max – ufl-compatible, the maximum velocity
dt – ufl-compatible, the time step

See:: https://en.wikipedia.org/wiki/Courant%E2%80%93Friedrichs%E2%80%93Lewy_condition

property timescheme

property t

property dt

set_initial_condition(u0, v0, t0=0) → None[source]

Apply initial conditions

Parameters:

u0 – u at t0
v0 – du/dt at t0
t0 – start time (default: 0)

u0 and v0 can be:

Python callable -> will be evaluated at nodes
-> e.g. u0 = lambda x: np.zeros((dim, x.shape[1]), dtype=PETSc.ScalarType)
scalar (int, float, complex, PETSc.ScalarType)
-> e.g. u0 = 0
array (list, tuple, np.ndarray) or fem.function.Constant
-> e.g. u0 = [0,0,0]
fem.function.Function
-> e.g. u0 = fem.Function(V)

solve(num_steps, **kwargs)[source]

class elastodynamicsx.solvers.EigenmodesSolver(comm: Comm, M: Mat, C: Mat, K: Mat, **kwargs)[source]

Bases: EPS

Convenience class inhereted from SLEPc.EPS, with methods and default parameters that are relevant for computing the resonances of an elastic component.

Parameters:

comm – The MPI communicator
M – The mass matrix
C – The damping matrix. C=None means no dissipation. C!=None is not supported yet.
K – The stiffness matrix

Keyword Arguments:

nev – The number of eigenvalues to be computed

Example

# ###
# Free resonances of an elastic cube
# ###
from mpi4py import MPI
from dolfinx import mesh, fem
from elastodynamicsx.solvers import EigenmodesSolver
from elastodynamicsx.pde import material, PDE

domain = mesh.create_box(MPI.COMM_WORLD, [[0,0,0], [1,1,1]], [10,10,10])
V      = dolfinx.fem.functionspace(domain, ("CG", 1, (3,)))

rho, lambda_, mu = 1, 2, 1
mat = material(V, rho, lambda_, mu)
pde = PDE(V, materials=[mat])
nev = 6 + 6  # the first 6 resonances are rigid body motion
eps = EigenmodesSolver(V.mesh.comm, pde.M(), None, pde.K(), nev=nev)
eps.solve()
eps.plot(V)
freqs = eps.getEigenfrequencies()
print('First resonance frequencies:', freqs)

getWn() → ndarray[source]: The eigen angular frequencies from the computed eigenvalues

getEigenfrequencies() → ndarray[source]: The eigenfrequencies from the computed eigenvalues

getEigenmodes(which='all') → List[Vec][source]

Returns the desired modeshapes

Parameters:: which – ‘all’, or an integer, or a list of integers, or a slice object

Example

getEigenmodes()   # returns all computed eigenmodes
getEigenmodes(3)  # returns mode number 4
getEigenmodes([3,5])  # returns modes number 4 and 6
getEigenmodes(slice(0,None,2))  # returns even modes

getModalBasis() → ModalBasis[source]

getErrors() → ndarray[source]: Returns the error estimate on the computed eigenvalues

plot(function_space: FunctionSpace, which='all', **kwargs) → None[source]

Plots the desired modeshapes

Parameters:

function_space – The underlying function space
which – ‘all’, or an integer, or a list of integers, or a slice object -> the same as for getEigenmodes

printEigenvalues() → None[source]: Prints the computed eigenvalues and error estimates

class elastodynamicsx.solvers.FrequencyDomainSolver(comm: Comm, M: Mat, C: Mat, K: Mat, b: Vec, b_update_function: Callable | None = None, **kwargs)[source]

Bases: object

Class for solving frequency domain problems.

Parameters:

comm – The MPI communicator
M – The mass matrix
C – The damping matrix
K – The stiffness matrix
b – The load vector
b_update_function – A function that updates the load vector (in-place) The function must take b,omega as parameters. e.g.: b_update_function = lambda b,omega: b[:]=omega If set to None, the call is ignored.

Keyword Arguments:

petsc_options – Options that are passed to the linear algebra backend PETSc. For available choices for the ‘petsc_options’ kwarg, see the PETSc documentation.

Example

from mpi4py import MPI

from dolfinx import mesh, fem
import ufl

from elastodynamicsx.solvers import FrequencyDomainSolver
from elastodynamicsx.pde import material, BoundaryCondition, PDE

# domain
length, height = 10, 10
Nx, Ny = 10, 10
domain = mesh.create_rectangle(MPI.COMM_WORLD, [[0,0], [length,height]], [Nx,Ny])
V      = dolfinx.fem.functionspace(domain, ("Lagrange", 1, (2,)))

# material
rho, lambda_, mu = 1, 2, 1
mat = material(V, rho, lambda_, mu)

# absorbing boundary condition
Z_N, Z_T = mat.Z_N, mat.Z_T  # P and S mechanical impedances
bcs = [ BoundaryCondition(V, 'Dashpot', (Z_N, Z_T)) ]

# gaussian source term
F0     = fem.Constant(domain, PETSc.ScalarType([1, 0]))  # polarization
R0     = 0.1  # radius
x0, y0 = length/2, height/2  # center
x      = ufl.SpatialCoordinate(domain)
gaussianBF = F0 * ufl.exp(-((x[0]-x0)**2 + (x[1]-y0)**2) /2 /R0**2) / (2 * 3.141596*R0**2)
bf         = BodyForce(V, gaussianBF)

# PDE
pde = PDE(V, materials=[mat], bodyforces=[bf], bcs=bcs)

# solve
fdsolver = FrequencyDomainSolver(V.mesh.comm, pde.M(), pde.C(), pde.K(), pde.b())
omega    = 1.0
u        = fem.Function(V, name='solution')
fdsolver.solve(omega=omega, out=u.x.petsc_vec)

default_petsc_options = {'ksp_type': 'preonly', 'pc_type': 'lu'}

solve(omega: float | ndarray, out: Vec | None = None, callbacks: List[Callable] = [], **kwargs) → Vec[source]

Solve the linear problem

Parameters:

omega – The angular frequency (scalar or array)
out – The solution (displacement field) to the last solve. If None a new PETSc.Vec is created
callbacks – If omega is an array, list of callback functions to be called after each solve (e.g. plot, store solution, …). Ignored if omega is a scalar.

Keyword Arguments: