Elastodynamics

This tutorial adds time dependence to the static elasticity problem. It introduces displacement and velocity as state fields and advances them with the public IDAAdapter execution path. IDA solves the canonical residual F(t, y, ydot) = 0; the application accepts the resulting state before native output is written.

The executable is elastodynamics, from apps/app_elastodynamics.cc. The canonical 2D input is tutorials/elastodynamics/strong_dirichlet.prm.in:

# First-order standalone adaptation of tutorials/elasticity/strong_dirichlet.prm.
set dimension = 2
set space dimension = 2

subsection Elastodynamics
  set FE degree                  = 1
  set Output directory            = @TEST_OUTPUT_DIR@/tutorial-output/elastodynamics-strong
  set Output name                 = strong_dirichlet
  set Number of refinement cycles = 5
  set Dirichlet boundary ids      = 0,1,2,3

  subsection Domain
    set Initial refinement          = 1
    set Grid generator           = hyper_cube
    set Grid generator arguments = 0: 1: true
    set Triangulation type       = distributed
  end

  subsection Material properties
    subsection default
      set Density       = 1.0
      set Lame mu       = 1.0
      set Lame lambda   = 2.0
      set Rayleigh alpha = 0
      set Rayleigh beta = 0
      set Viscosity eta = 0
    end
  end

  subsection Time interval
    set Initial time       = 0.0
    set Final time         = 0.01
    set Output time interval = 0.01
  end
  subsection Fixed step
    set Time step           = 0.01
    set Number of time steps = 1
    set Policy               = number_of_steps
  end

  subsection Functions
    subsection Body force
      set Function constants  =
      set Function expression = -2*x^2 - 12*x*y + 8*x - 8*y^2 + 14*y + 5*pi^2*sin(pi*x)*sin(pi*y) - 3; -8*x^2 - 12*x*y + 14*x - 2*y^2 + 8*y - 3*pi^2*cos(pi*x)*cos(pi*y) - 3
      set Variable names      = x,y,t
    end
    subsection Displacement boundary
      set Function constants  =
      set Function expression = x*y*(x - 1)*(y - 1) + sin(pi*x)*sin(pi*y); x*y*(x - 1)*(y - 1)
      set Variable names      = x,y,t
    end
    subsection Velocity boundary
      set Function constants  =
      set Function expression = 0; 0
      set Variable names      = x,y,t
    end
    subsection Initial displacement
      set Function constants  =
      set Function expression = x*y*(x - 1)*(y - 1) + sin(pi*x)*sin(pi*y); x*y*(x - 1)*(y - 1)
      set Variable names      = x,y,t
    end
    subsection Initial velocity
      set Function constants  =
      set Function expression = 0; 0
      set Variable names      = x,y,t
    end
    subsection Exact solution
      set Function constants  =
      set Function expression = x*y*(x - 1)*(y - 1) + sin(pi*x)*sin(pi*y); x*y*(x - 1)*(y - 1)
      set Variable names      = x,y,t
    end
  end

  subsection Solver
    subsection Control
      set Log frequency = 1
      set Log history   = false
      set Log result    = false
      set Max steps     = 500
      set Reduction     = 1.e-12
      set Tolerance     = 1.e-12
    end
  end
end

subsection Error
  set Enable computation of the errors = true
  set Error file name                  =
  set Error precision                  = 8
  set Exponent for p-norms             = 2
  set Extra columns                    = cells, dofs
  set List of error norms to compute   = L2_norm, Linfty_norm, H1_norm
  set Rate key                         = dofs
  set Rate mode                        = reduction_rate_log2
end

The application reads dimension and space dimension from the file. This input selects the full-dimensional 2D solver; the executable also supports the full-dimensional 3D instantiation. The input uses one time step and a manufactured displacement so it stays a small, deterministic smoke example.

Run it with:

cmake --build build -j
./build/elastodynamics_debug \
  build/tutorials/elastodynamics/strong_dirichlet.prm

Output is written below build/test_output/tutorial-output/elastodynamics-strong. The same input is exercised by AppExecutables.TutorialElastodynamics.

The physical model exposes separate mass, stiffness, and damping operators. The application registers the Problem directly with IDAAdapter, marks both displacement and velocity as differential Fields, and installs an accepted output callback. initial_acceleration() supplies a physically consistent initial derivative; accept_state() updates the Problem only after IDA has accepted or interpolated a state. Solver-neutral residual and execution-adapter concepts are described in Architecture concepts. The larger convergence cases remain in tutorials/elastodynamics/ and are listed in the verification reference.

Materials are selected under Elastodynamics/Material properties. Map mesh material ids to named subsections with Material tags by material id, for example 0: default, 1: stiff. Unmapped ids use default; several ids may share a tag. Each material accepts Density, Lame mu, Lame lambda, Rayleigh alpha, Rayleigh beta, and Viscosity eta.

The damping bilinear form on each cell is

\[ C(v,w) = \alpha\rho (v,w) + \beta\bigl(2\mu (\varepsilon(v),\varepsilon(w)) + \lambda (\nabla\cdot v,\nabla\cdot w)\bigr) + \eta (\varepsilon(v),\varepsilon(w)). \]

Density and shear modulus must be finite and positive. The remaining coefficients must be finite and non-negative. The viscosity convention is the same as in ElasticityProblem: there is no factor of two before Viscosity eta. To migrate an older elastodynamics input, move the Material entries to Material properties/default and replace Damping shear = s with Viscosity eta = 2*s. The former Damping bulk parameter has no independent counterpart in this schema; Rayleigh beta scales the entire elastic stiffness.

ElastodynamicsSolver also runs on simplex meshes (triangles in 2D, tetrahedra in 3D). Read the mesh from a file and select the fullydistributed triangulation, since p4est holds only quadrilaterals and hexahedra:

subsection Elastodynamics
  subsection Domain
    set Grid generator = liver.msh
    set Triangulation type = fullydistributed
  end
end

On simplex meshes the solver uses FE_SimplexP, QGaussSimplex, and a linear MappingFE; meshes that mix simplex and hypercube cells are rejected. Initial refinement refines the mesh before it is distributed. The fullydistributed triangulation cannot be refined afterwards, so Number of refinement cycles must stay at 1. The ParsedConvergenceTable error integration supports only hypercube cells: run() skips the error table on simplex meshes and compute_error() throws. gtests/elastodynamics_simplex_01.cc runs a 48-tetrahedron cube from data/tests/simplex_cube.msh.in with nonzero body force and Dirichlet data, in serial and on two MPI ranks. The generated gtests/parameters/elastodynamics_simplex.prm also checks quadratic spatial displacement and cubic time dependence with the native trapezoidal method.

The elastodynamics_adapter executable accepts the same uniform simplex meshes through Elastodynamics/Domain. Set Triangulation type = fullydistributed and use one refinement cycle. Its IDA solve uses simplex operators and skips the convergence table. The generated test input gtests/parameters/elastodynamics_adapter_simplex.prm supplies nonzero body force and time-dependent displacement and velocity boundary values.

Next, Reduced Poisson introduces a lower-dimensional geometry and coupling.

Robin tractions

A displacement spring on the boundary satisfies

\[ \sigma(u,v)n + \beta_R(x,t)u = g(x,t), \qquad \beta_R(x,t)\geq 0, \]

where \(v=\dot u\) and \(\sigma(u,v)\) includes the elastic and viscous stresses of the material model above. The scalar Robin coefficient \(\beta_R\) is independent of the material’s Rayleigh beta. It multiplies displacement; it does not introduce boundary damping proportional to velocity.

For a colored unit square, this example prescribes displacement on three faces and a time-dependent Robin traction on the right face:

subsection Elastodynamics
  set Dirichlet boundary ids = 0, 2, 3
  set Neumann boundary ids =
  set Robin boundary ids = 1
  subsection Functions
    subsection Robin coefficient
      set Function expression = 2+x+t
    end
    subsection Robin boundary
      set Function expression = sin(t); 0
    end
  end
end

Set Domain/Grid generator arguments = 0: 1: true for that face numbering. The vector load has one component per ambient spatial coordinate. Coefficients must remain finite and nonnegative, loads must be finite, and Robin ids must be disjoint from displacement and Neumann ids.

Both the native fixed-step solver and the IDA adapter evaluate the coefficient and load at the requested physical time. In the native trapezoidal method, the previous-time elastic force uses the previous coefficient. Time-varying loads can be combined with a constant boundary spring. A Robin update leaves the mass and material damping operators unchanged.

The same entries are available in the matrix and fiber Elastodynamics subsections of fiber_reinforced_elastodynamics, and in the solid subsection of metric_flow_x_elastodynamics. Hypercube and simplex face integration use the same reference-cell choice as the volume operators.