Harmonic reinforcement with or without a lift

The matrix always has dimension spacedim. The root dimension parameter is now the dimension of the reinforcement mesh; space dimension is the ambient dimension. The application accepts all six pairs with 1 <= dimension <= space dimension <= 3.

• If dim < spacedim, the reinforcement is lifted through a cross section of dimension spacedim-dim to obtain a volume.

• If dim == spacedim, its independently supplied volume mesh is used directly. There is no Lift subsection in this case.

Fiber Reinforced Elasticity Frequency/Formulation selects the displacement model. Both choices use the reinforcement’s excess material coefficients, which are added to the matrix response already present in that volume.

Run the examples

Configure with MPI’s C++ compiler, then build the application:

cmake -S . -B build-debug -DCMAKE_BUILD_TYPE=Debug \
  -DCMAKE_CXX_COMPILER=mpicxx -DDEAL_II_DIR=/path/to/deal.II \
  -DENABLE_COUPLED_PROBLEMS=OFF -DIMMERSX_BUILD_CORAL_PLUGINS=OFF
cmake --build build-debug --target fiber_reinforced_elasticity_freq_debug -j 4
mpirun -np 2 build-debug/fiber_reinforced_elasticity_freq_debug \
  build-debug/tutorials/fiber_reinforced_elasticity_freq/parameters.prm

The Release target is fiber_reinforced_elasticity_freq. CMake configures the parameter templates into the build tree, with absolute output paths. The executable and input can also be addressed by absolute paths from another working directory.

Six complete inputs are provided under tutorials/fiber_reinforced_elasticity_freq/:

Generated input	Reinforcement	Formulation
parameters.prm	1D in 2D, with a 1D lift	independent
shared.prm	1D in 2D, with a 1D lift	shared
full_dimension.prm	2D in 2D, without a lift	independent
fiber_3d.prm	1D in 3D, with a 2D lift	independent
fiber_3d_shared.prm	1D in 3D, with a 2D lift	shared
sheet.prm	2D in 3D, with a 1D lift	independent

The fiber in a 3D matrix tutorial explains the codimension-two geometry and both displacement formulations.

The canonical lifted input is:

# Matrix and additive/excess reinforcement under harmonic traction.
set dimension = 1
set space dimension = 2

subsection Fiber Reinforced Elasticity Frequency
  set Formulation = independent
  set Frequencies = 0, 0.1, 0.2
  set Output directory = @TEST_OUTPUT_DIR@/tutorial-output/fiber-reinforced-elasticity-frequency
  set Multiplier FE degree = 0

  subsection Lift
    set Thickness = 0.1
  end

  subsection Solver
    set Solver = direct
  end

  subsection Matrix
    set FE degree = 1
    set Initial refinement = 2
    set Dirichlet boundary ids = 0
    set Neumann boundary ids = 1
    subsection Grid generation
      set Grid generator = hyper_cube
      set Grid generator arguments = -1: 1: true
      set Triangulation type = distributed
    end
    subsection Material properties
      subsection default
        set Lame mu = 1
        set Lame lambda = 1
        set Density = 1
        set Viscosity eta = 0.05
        set Rayleigh alpha = 0.1
        set Rayleigh beta = 0.001
      end
    end
    subsection Frequency domain
      set Neumann phase = 15
    end
    subsection Functions
      subsection Right hand side
        set Function expression = 0; 0
      end
      subsection Dirichlet boundary conditions
        set Function expression = 0; 0
      end
      subsection Neumann boundary conditions
        set Function expression = 1; 0
      end
    end
  end

  subsection Fiber
    set FE degree = 1
    set Initial refinement = 1
    set Dirichlet boundary ids =
    set Neumann boundary ids =
    subsection Grid generation
      set Grid generator = hyper_cube
      set Grid generator arguments = -0.65: 0.65: true
      set Triangulation type = fullydistributed
    end
    subsection Material properties
      subsection default
        # Additive/excess coefficients per unit reinforcement volume.
        set Lame mu = 8
        set Lame lambda = 8
        set Density = 0.2
        set Viscosity eta = 0.05
        set Rayleigh alpha = 0.1
        set Rayleigh beta = 0.001
      end
    end
    subsection Functions
      subsection Right hand side
        set Function expression = 0; 0
      end
    end
  end
end

Each example clamps the matrix’s left side, applies a traction with a phase of 15 degrees on the right, and sweeps 0, 0.1, and 0.2 Hz with nonzero damping. Change Formulation to switch models. Use a separate output directory for each run so old datasets from another model are not mistaken for current results.

Iterative solution

For shared displacement, select the positive elasticity surrogate with:

subsection Fiber Reinforced Elasticity Frequency
  set Formulation = shared
  subsection Solver
    set Solver = iterative
    set Preconditioner = block_diagonal
  end
end

The surrogate is \(P_0=K+\omega^2M+\omega D\). Its AMG hierarchy uses elastic rigid-body modes and is shared by the real and imaginary displacement fields. Shared-displacement solves include both the matrix and the retained reinforcement material in \(P_0\). Each application uses at most 30 AMG-preconditioned CG iterations with a relative residual target of 0.1. The outer FGMRES checks the residual of the original harmonic system against Solver/Tolerance.

Preconditioner = none remains available. A positive surrogate does not remove physical resonances; weakly damped cases may need more outer iterations.

Composed normal preconditioner

normal currently requires the shared formulation. Select it with:

subsection Fiber Reinforced Elasticity Frequency
  set Formulation = shared
  subsection Solver
    set Solver = iterative
    set Preconditioner = normal
    set Tolerance = 1e-7
    set FGMRES basis size = 100
    set Normal inner maximum iterations = 100
    set Normal inner reduction = 1e-3
  end
end

Writing \(\mathcal H\) for the assembled real/imaginary system and \(\mathcal P=\operatorname{diag}(P_0,P_0)\), this applies \(\mathcal P^{-1}\mathcal H^T\mathcal P^{-1}\) as a right preconditioner. Each application uses four real CG/AMG solves and one transpose matrix action. The constrained rows use the actual diagonal of the corresponding harmonic block. The physical matrix and right-hand side are unchanged.

The inner CG stops at its relative reduction target or iteration cap. FGMRES checks the residual of the original system, including a final true-residual check. High material contrast can require more work in the inner solves; Normal inner reduction and Normal inner maximum iterations control that tradeoff. These two settings apply to normal; the existing block_diagonal policy retains its reduction of 0.1 and cap of 30.

MUMPS-BLR preconditioner

With a MUMPS-enabled deal.II build and the Trilinos backend, use an approximate factorization of the complete system inside FGMRES:

subsection Fiber Reinforced Elasticity Frequency
  subsection Solver
    set Solver = iterative
    set Preconditioner = mumps_blr
    set Tolerance = 1e-7
    set FGMRES basis size = 100
    subsection MUMPS BLR
      set Tolerance = 1e-5
      set Memory limit MB per rank = 0
      set Output details = true
    end
  end
end

MUMPS BLR/Tolerance is the absolute truncation threshold after scaling frontal matrices; it is separate from the outer Tolerance. A larger threshold reduces the factor storage but can weaken the preconditioner. The factorization uses nonsymmetric real LU, MUMPS scaling 8 and UCFS; FGMRES supplies the final accuracy instead of internal MUMPS refinement. The factors are reused for every Krylov application within the solve.

The memory setting is a MUMPS workspace budget in MB per rank, not a limit on the complete process RSS. Zero leaves that budget unrestricted. Output details enables MUMPS analysis, factorization and memory diagnostics. The outer solve rejects invalid corrections and checks the true residual. The current MUMPS vector interface requires the global system size to be at least the MPI rank count; the adapter redistributes its output to the semantic field partitions, including ranks without original field entries.

The same LinearAdapter policies can be used by custom contributors. normal requires a registered local inverse for each field; factories receive that field’s current diagonal matrix and vector reinitializer. It composes those inverses with dealii::transpose_operator and operator products. mumps_blr requires an assemblable global matrix and does not require local inverses.

Lift geometry and material units

The Lift subsection is registered only for reduced reinforcement meshes. Its cross section has dimension spacedim-dim. The supported hyper_cube section is [-1,1]^(spacedim-dim), scaled by Thickness = h; its physical measure is therefore (2h)^(spacedim-dim). This is an interval in codimension one and a square in codimension two. A 2D sheet in 3D extrudes along its local normal, using both tangent directions.

Lift/Representative quadrature and Lift/Cross section control the quadrature points and repetitions. Increase them to resolve the background FE traces and material interfaces. Points are located in the distributed matrix mesh through deal.II’s particle facilities; every lifted point must lie inside the matrix’s searchable domain. Quadrature is not split at matrix-cell interfaces, so increasing only the Gauss order need not give monotone convergence.

Material coefficients are per unit physical volume. The lift supplies the section measure; do not multiply the material coefficients by it again. The current harmonic reconstruction uses the constant section mode: Maximum basis degree = 0 and Selected modes = 0. The thickness is a positive numeric constant. The reduced harmonic model assumes flat cell manifolds; curved geometries can be approximated with piecewise-flat meshes. Geometry remains fixed during a sweep. The quadrature uses representative measure times section measure; curvature corrections for offset Jacobians and derivatives of transported frames are not included.

Frequency and phase convention

For \(u(t)=\operatorname{Re}(\widehat u\exp(i\omega t))\), each material contributes

\[ H(\omega)=K-\omega^2M+i\omega D,\qquad \omega=2\pi f, \]

with

\[\begin{split} \begin{aligned} K(u,v)&=\int 2\mu\,\varepsilon(u):\varepsilon(v) +\lambda\,\operatorname{div}u\operatorname{div}v,\\ M(u,v)&=\int\rho\,u\cdot v,\\ D(u,v)&=\beta K(u,v)+\alpha M(u,v) +\int\eta\,\varepsilon(u):\varepsilon(v). \end{aligned} \end{split}\]

The real matrix representation is

\[\begin{split} \begin{pmatrix}K-\omega^2M&-\omega D\\ \omega D&K-\omega^2M\end{pmatrix}. \end{split}\]

The root Frequencies list is shared by matrix and reinforcement. Their Material properties and Functions subsections remain independent. Spatial load amplitudes are evaluated at t=0, then combined with phase angles in degrees. Reinforcement material ids are taken from its own mesh, including in the shared formulation.

Independent displacements

Formulation = independent retains separate matrix and reinforcement Fields, with real and imaginary multiplier constraints. For reduced reinforcement, the lift reconstructs each ambient displacement component as constant across the section. Its native tangential strain, mass, damping, body forces, and endpoint/boundary tractions are weighted by the physical section measure. A full-dimensional reinforcement uses its ordinary volumetric strain and volume measure.

Writing \(E\) for constant normal extension, the continuity equations are

\[ \int_{\Omega_f}E\lambda\cdot(\widehat u_m-E\widehat u_f)=0. \]

This is a weak section-average constraint. It does not require every matrix point in a section to have the same displacement. For full-dimensional reinforcement, \(E\) is the identity. Both real and imaginary equations are imposed. The complex system is

\[\begin{split} \begin{pmatrix} H_m&0&C_m^T\\ 0&H_f&-C_f^T\\ C_m&-C_f&0 \end{pmatrix} \begin{pmatrix}\widehat u_m\\\widehat u_f\\\widehat\lambda\end{pmatrix} = \begin{pmatrix}\widehat f_m\\\widehat f_f\\0\end{pmatrix}. \end{split}\]

The multiplier defaults to vector-valued continuous FE_Q on the reinforcement mesh. Multiplier FE degree overrides its degree; Multiplier space/Finite element selects an explicit element with spacedim components. The mesh and space must yield independent constraints. The default direct solver treats the complete saddle-point system without assuming positive definite or individually invertible displacement blocks. An undamped resonance may still make the complete system singular.

This is a vector small-strain model with constant section reconstruction. It has no independent rotations or higher section modes.

Shared displacement

Formulation = shared registers only the matrix’s real and imaginary spatial displacement Fields. The reinforcement mesh supplies integration geometry and material ids; it carries no independent displacement unknowns and introduces no multiplier. The equation is

\[ (H_m+H_f^{\mathrm{bulk}}+B_{\beta,f})\widehat u =\widehat f_m+\widehat f_f^{\mathrm{bulk}}+\widehat g_f. \]

Reinforcement stiffness, mass, and damping use the matrix displacement’s full ambient gradients at physical reinforcement quadrature points. Its body force is evaluated at those same points, with its own phase and material overrides. Quadrature points outside the matrix domain are omitted from reinforcement stiffness, mass, damping, and forcing. This selection uses the physical lifted points, so a partially overlapping cross section retains its interior samples. A completely exterior reinforcement contributes no material or load.

Fiber/Robin boundary ids selects springs and applied tractions on the reinforcement boundary. For a reduced reinforcement, each selected representative face is extruded through the physical section: a rod endpoint gives an endcap, and a sheet edge gives an extruded edge surface. Full-dimensional reinforcement uses its own mesh faces. The spring term is

\[ v^T B_{\beta,f}u=\int_{\Gamma_f}\beta_f\,u\cdot v\,dS, \qquad v^T\widehat g_f=\int_{\Gamma_f}\widehat g_f\cdot v\,dS. \]

Boundary samples outside the matrix domain are omitted, using the same physical point selection as the reinforcement volume. Coefficients and loads are evaluated at the retained physical surface points. The scalar coefficient is finite, nonnegative and time independent; the traction has spacedim components. Their subsections are Fiber/Functions/Robin coefficient and Fiber/Functions/Robin boundary conditions. A numbered suffix selects an override for one boundary id. Fiber/Frequency domain/Robin phase and Robin phase by boundary id specify load phases in degrees. The spring acts on both real and imaginary displacement; its coefficient has no load phase.

In the independent formulation, the same entries apply to the independent fiber displacement. Reduced boundary data are evaluated on the representative mesh and weighted by the physical section measure, consistently with its constant normal reconstruction. Matrix Robin data use the corresponding subsections below Matrix in either formulation.

Essential, normal-flux and Neumann Fiber boundary-id sets remain empty for shared displacement; those constraints and tractions belong to the matrix. The default reinforcement boundary-id sets are empty.

This formulation allows the matrix displacement to vary across the lifted section. It is therefore a different reduced model from constant-extension independent displacements, even with the same coefficients and geometry.

Output and current limits

Both formulations write frequencies.csv and matrix VTU/PVTU datasets with real and imaginary displacement, component amplitudes, and phases in degrees. Independent displacements also write reinforcement datasets on its native mesh and multiplier datasets under interaction/frequency_<index>.

residual_norm() reports the complete algebraic residual. continuity_norm() reports the combined multiplier-row residual in the independent formulation and is zero when no constraint is present. reinforcement_measure() reports the integrated physical reinforcement volume. For shared displacement it includes only quadrature samples inside the matrix. These diagnostics do not estimate the discretization error.

Each Problem uses one prepared mesh during the frequency sweep. Set Number of refinement cycles = 1 and choose initial refinements before running. The adapter path supports strong Dirichlet conditions. Weak Dirichlet conditions remain available in standalone elasticity_freq; normal-flux constraints require a full-dimensional 2D or 3D Problem.

Compose a harmonic Problem

An uncoupled harmonic Problem can also be registered directly with LinearAdapter:

using namespace ImmersX;
ElasticityFreqProblemParameters<2> parameters;
initialize_parameters(parameter_file);
ElasticityFreqProblem<2> problem(parameters);
problem.make_grid();
problem.setup_fe();
problem.setup_dofs();
problem.prepare_frequency(1.0);

LinearAdapterParameters solver_parameters;
solver_parameters.solver = LinearSolver::direct;
LinearAdapter<ImmersXLA::MPI::Vector, ImmersXLA::MPI::BlockVector>
  adapter(solver_parameters, problem.communicator());
const auto fields = adapter.add(problem, "solid").fields();
auto state = adapter.make_state();
adapter.solve(state);
problem.accept_frequency_solution(adapter.field(state, fields.real),
                                  adapter.field(state, fields.imaginary));
problem.output_results(0);

Include immersx/physics/elasticity_freq_contributor.h and immersx/core/linear_adapter.h. The two algebraic Fields can be used in ordinary FE expressions and continuity constraints. A new Interaction can contribute to either residual row or introduce its own multiplier Field. Mesh ownership, material lookup, and essential constraints remain with the Problem; execution block numbers remain inside the adapter.

Each frequency has constant coefficients during one linear solve. Prepare the new frequency and create a fresh adapter composition for the next solve so its weak-term assembly actions use the new coefficients. For a new constitutive term, add its real and imaginary weak contributions with this same phase convention; frequency-dependent solver policy belongs in the execution adapter.