Overview

Introduction

ElastodynamiCSx deals with the following PDE:

\[\mathbf{M}\mathbf{a} + \mathbf{C}\mathbf{v} + \mathbf{K}(\mathbf{u}) = \mathbf{b},\]

where \(\mathbf{u}\), \(\mathbf{v}=\partial_t \mathbf{u}\), \(\mathbf{a}=\partial_t^2\mathbf{u}\) are the displacement, velocity and acceleration fields, \(\mathbf{M}\), \(\mathbf{C}\) and \(\mathbf{K}\) are the mass, damping and stiffness forms, and \(\mathbf{b}\) is an applied force. For time domain problems \(\mathbf{K}\) may be a non-linear function of \(\mathbf{u}\).

The library provides a high level interface to build and solve common problems in a few lines code.

Build problems

Using the elastodynamicsx.pde package:

• Common material laws:

  # V is a dolfinx.fem.FunctionSpace
  # cell_tags is a dolfinx.mesh.MeshTags object
  
  from elastodynamicsx.pde import material
  
  tag_mat1 = 1        # suppose this refers to cells associated with material no 1
  tags_mat2 = (2, 3)  # same for material no 2
  mat1 = material((V, cell_tags, tag_mat1), 'isotropic', rho=1, lambda_=2, mu=1)
  mat2 = material((V, cell_tags, tags_mat2), 'isotropic', rho=2, lambda_=4, mu=2)
  

  • linear elasticity: scalar, isotropic, cubic, hexagonal, trigonal, tetragonal, orthotropic, monoclinic, triclinic

    • damping laws: Rayleigh damping

  • hyperelastic: Saint Venant-Kirchhoff, Murnaghan

  • perfectly matched layers (PML): in the near future…

• Common boundary conditions (BCs):

  # facet_tags is a dolfinx.mesh.MeshTags object
  
  from elastodynamicsx.pde import boundarycondition
  
  tag_top = 1  # top boundary
  T_N  = fem.Constant(V.mesh, default_scalar_type([0, 1]))  # boundary load
  bc1  = boundarycondition((V, facet_tags, tag_top) , 'Neumann', T_N)  # prescribed traction
  

  • BCs involving \(\mathbf{u}\) and \(\boldsymbol{\sigma} . \mathbf{n}\): Free, Clamp, Dirichlet, Neumann, Robin

  • BCs involving \(\mathbf{v}\) and \(\boldsymbol{\sigma} . \mathbf{n}\): Dashpot

  • Multi-point constraint BCs: Periodic

• User-defined material laws and BCs, using the `ufl` library:

  Custom material laws

  Specify \(\mathbf{M}\), \(\mathbf{C}\) and \(\mathbf{K}\):

  import ufl
  
  # ###
  # Here we re-implement mat1 using the interface for custom material laws
  dx_mat1 = ufl.Measure("dx", domain=V.mesh, subdomain_data=cell_tags)(tag_mat1)
  
  # mass form
  m = lambda u, v: 1 * ufl.inner(u, v) * dx_mat1
  
  # stiffness form
  epsilon = lambda u: ufl.sym(ufl.grad(u))
  sigma = lambda u: 2 * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * 1 * epsilon(u)
  k = lambda u, v: ufl.inner(sigma(u), epsilon(v)) * dx_mat1
  
  mat1_user_defined = material(V, 'custom', is_linear=True, M_fn=m, K_fn=k)
  
  Custom BCs

  Specify \(\mathbf{C}\), \(\mathbf{K}\) and \(\mathbf{b}\):

  import ufl
  
  # ###
  # Here we re-implement bc1 using the interface for custom BCs
  ds_bc1 = ufl.Measure("ds", domain=V.mesh, subdomain_data=facet_tags)(tag_top)
  
  # right hand side term
  b = lambda v: ufl.inner(T_N, v) * ds_bc1
  
  bc1_user_defined = boundarycondition(V , 'custom', b_fn=b)
  

• Define body forces:

  from elastodynamicsx.pde import BodyForce
  
  amplitude = default_scalar_type([0, 0])  # a dummy load amplitude
  def shape_x(x):
      x1, x2 = 0.2, 0.3
      y1, y2 = 0.4, 0.5
      return (x[0] >= x1) * (x[0] <= x2) * (x[1] >= y1) * (x[1] <= y2)  # a dummy shape
  
  f_body = fem.Function(V)
  f_body.interpolate(lambda x: amplitude[:, np.newaxis] * shape_x(x)[np.newaxis, :])
  f1 = BodyForce((V, cell_tags, None), f_body)  # None for the entire domain
  

• Assemble several materials, BCs and body forces into a PDE instance:

  from elastodynamicsx.pde import PDE
  
  pde = PDE(V, materials=[mat1, mat2], bodyforces=[f1], bcs=[bc1])
  
  # M, C, K, b form functions: pde.M_fn, pde.C_fn, pde.K_fn, pde.b_fn
  # eigs / freq. domain -> M, C, K matrices:    pde.M(),  pde.C(),  pde.K()
  # waveguides          -> K0, K1, K2 matrices: pde.K0(), pde.K1(), pde.K2()
  

  • Get the \(\mathbf{M}\), \(\mathbf{C}\), \(\mathbf{K}\) weak forms - ufl format

  • Compile the \(\mathbf{M}\), \(\mathbf{C}\), \(\mathbf{K}\) matrices - petsc format

• Build the weak form of a time domain problem

  • Explicit schemes: leapfrog

  • Implicit schemes (currently restricted to linear PDEs): Newmark-beta, midpoint, linear acceleration, HHT-alpha, generalized-alpha

Solve problems

Using the elastodynamicsx.solvers package:

Time domainFrequency domainEigenmodesGuided waves

# Time integration
from elastodynamicsx.solvers import TimeStepper

dt, num_steps = 0.01, 100  # t=[0..1)

# Define a function that will update the source term at each time step
def update_T_N_function(t):
    forceVector = default_scalar_type([0, 1])
    T_N.value   = np.sin(t) * forceVector

# Initialize the time stepper: compile forms and assemble the mass matrix
tStepper = TimeStepper.build(V,
                             pde.M_fn, pde.C_fn, pde.K_fn, pde.b_fn, dt, bcs=pde.bcs,
                             scheme='newmark')

# Define the initial values
tStepper.set_initial_condition(u0=[0, 0], v0=[0, 0], t0=0)

# Solve: run the loop on time steps; live-plot the result every 10 steps
tStepper.solve(num_steps-1,
               callfirsts=[update_T_N_function],
               callbacks=[],
               live_plotter={'refresh_step':10, 'clim':[-1, 1]})

Post-process solutions

Using the elastodynamicsx.solutions package:

• Eigenmodes solutions:

# eps is a elastodynamicsx.solvers.EigenmodesSolver
# eps.solve() has already been performed

# Get the solutions
mbasis = eps.getModalBasis()  # a elastodynamicsx.solutions.ModalBasis

# Access data
eigenfreqs = mbasis.fn     # a np.ndarray
modeshape5 = mbasis.un[5]  # a PETSc.Vec vector

# Visualize
mbasis.plot(function_space=V)  # V is a dolfinx.fem.FunctionSpace