Variational form and the discretised equations
==============================================

This chapter develops the linear-elastic finite-element
equations femorph-solver assembles and solves.  The treatment
follows Bathe (2014) §4 and Zienkiewicz & Taylor (2013) §2 +
§8 closely; readers familiar with FE foundations can skip to
:doc:`isoparametric` for the geometric mapping.

Strong form: Cauchy momentum
----------------------------

For a linear-elastic body :math:`\Omega \subset \mathbb{R}^{3}`
with displacement field :math:`\mathbf{u}(\mathbf{x})`, body
force :math:`\mathbf{b}(\mathbf{x})`, and surface traction
:math:`\mathbf{t}(\mathbf{x})` prescribed on
:math:`\partial \Omega_{t} \subseteq \partial \Omega`, static
equilibrium reads

.. math::
    :label: strong-form

    \nabla \cdot \boldsymbol{\sigma} + \mathbf{b} = \mathbf{0}
    \quad \text{in } \Omega,

with the constitutive law :math:`\boldsymbol{\sigma} = \mathbf{C}
\boldsymbol{\varepsilon}` (Hooke) and the small-strain kinematics

.. math::

    \boldsymbol{\varepsilon} = \tfrac{1}{2}\!
    \bigl(\nabla \mathbf{u} + (\nabla \mathbf{u})^{\!\top}\bigr).

Boundary conditions partition :math:`\partial \Omega` into

* a Dirichlet portion :math:`\partial \Omega_{u}` where
  :math:`\mathbf{u} = \bar{\mathbf{u}}` is prescribed, and
* a Neumann portion :math:`\partial \Omega_{t}` where
  :math:`\boldsymbol{\sigma}\,\hat{\mathbf{n}} = \mathbf{t}`
  is prescribed,

with :math:`\partial \Omega_{u} \cap \partial \Omega_{t} =
\varnothing`.  The dynamic generalisation adds the inertia
term :math:`\rho\, \ddot{\mathbf{u}}` on the left.

Weak form: principle of virtual work
------------------------------------

Multiplying :eq:`strong-form` by a test function
:math:`\delta \mathbf{u}` that vanishes on
:math:`\partial \Omega_{u}`, integrating over :math:`\Omega`,
and using integration by parts (Gauss–Ostrogradsky) on the
divergence term gives

.. math::
    :label: weak-form

    \int_{\Omega}
        \boldsymbol{\varepsilon}(\delta \mathbf{u})^{\!\top}
        \mathbf{C}
        \boldsymbol{\varepsilon}(\mathbf{u})
        \,\mathrm{d}V
    \;=\;
    \int_{\Omega}
        \delta \mathbf{u}^{\!\top}\, \mathbf{b}
        \,\mathrm{d}V
    +
    \int_{\partial \Omega_{t}}
        \delta \mathbf{u}^{\!\top}\, \mathbf{t}
        \,\mathrm{d}A.

Equation :eq:`weak-form` is the **principle of virtual work**:
the internal virtual strain energy (left) balances the external
virtual work (right) for every admissible
:math:`\delta \mathbf{u}`.  No second derivatives appear, so the
trial / test spaces only need :math:`H^{1}(\Omega)` regularity
— a key relaxation that admits piecewise-polynomial finite-
element approximations.

Galerkin discretisation
-----------------------

Approximate :math:`\mathbf{u}` and :math:`\delta \mathbf{u}` in
the **same** finite-dimensional space (Bubnov-Galerkin):

.. math::

    \mathbf{u}(\mathbf{x}) \approx
    \mathbf{N}(\mathbf{x})\, \mathbf{u}_e,
    \qquad
    \delta \mathbf{u}(\mathbf{x}) \approx
    \mathbf{N}(\mathbf{x})\, \delta \mathbf{u}_e,

where :math:`\mathbf{N}` is a per-element shape-function matrix
of size :math:`(\dim, n_\mathrm{dof,e})` and :math:`\mathbf{u}_e`
collects the element-local nodal DOFs.  The strain follows from
the strain-displacement matrix
:math:`\mathbf{B}(\mathbf{x})`,

.. math::

    \boldsymbol{\varepsilon}(\mathbf{u}) =
    \mathbf{B}(\mathbf{x})\, \mathbf{u}_e,

with :math:`\mathbf{B}` containing the spatial derivatives of
:math:`\mathbf{N}`.  Substituting into :eq:`weak-form` and
requiring the equation to hold for every
:math:`\delta \mathbf{u}_e` yields the per-element bilinear
forms

.. math::
    :label: ke-me

    \mathbf{K}_e =
    \int_{\Omega_e}
        \mathbf{B}^{\!\top} \mathbf{C}\, \mathbf{B}
        \,\mathrm{d}V,
    \qquad
    \mathbf{M}_e =
    \int_{\Omega_e}
        \rho\, \mathbf{N}^{\!\top}\, \mathbf{N}
        \,\mathrm{d}V,
    \qquad
    \mathbf{f}_e =
    \int_{\Omega_e}
        \mathbf{N}^{\!\top}\, \mathbf{b}
        \,\mathrm{d}V
    +
    \int_{\partial \Omega_e \cap \partial \Omega_t}
        \mathbf{N}^{\!\top}\, \mathbf{t}
        \,\mathrm{d}A.

Assembling per-element contributions into the global vectors
gives the linear-elastic static problem

.. math::
    :label: global-system

    \mathbf{K}\, \mathbf{U} = \mathbf{F},

and its dynamic counterpart

.. math::
    :label: dynamic-system

    \mathbf{M}\, \ddot{\mathbf{U}} +
    \mathbf{K}\, \mathbf{U} = \mathbf{F}(t),

from which the modal eigenvalue problem
:math:`\mathbf{K}\, \boldsymbol{\phi} = \omega^{2}\,
\mathbf{M}\, \boldsymbol{\phi}` falls out by separation of
variables.

How femorph-solver evaluates each integral
------------------------------------------

* :math:`\mathbf{B}` is supplied by the element kernel
  (``HEX8.ke``, ``TET10.ke``, …) — see
  :doc:`../elements/index`.
* :math:`\mathbf{C}` is built from the element's material via
  :func:`femorph_solver.elements._stress.elasticity_matrix` for
  isotropic linear elasticity (Lamé constants from
  :math:`E, \nu`).
* :math:`\rho` is read off the material's ``DENS`` field.
* :math:`\mathrm{d}V` is computed via the isoparametric
  mapping (:doc:`isoparametric`) and the integral is
  evaluated by Gauss quadrature on the reference element.
* :math:`\mathbf{f}_e` is currently assembled from explicit
  nodal forces (``Model.apply_force``) — distributed-pressure
  surface integrals land alongside the PLOAD4 / DSLOAD interop
  work (TA-11).

Reduction with Dirichlet BCs
----------------------------

Prescribed-displacement DOFs (`Model.fix(node, dof, value)`)
are eliminated by **row/column partitioning** of
:math:`\mathbf{K}` into free / constrained blocks:

.. math::

    \begin{bmatrix}
        \mathbf{K}_{ff} & \mathbf{K}_{fc} \\
        \mathbf{K}_{cf} & \mathbf{K}_{cc}
    \end{bmatrix}
    \begin{bmatrix}
        \mathbf{u}_f \\ \mathbf{u}_c
    \end{bmatrix}
    =
    \begin{bmatrix}
        \mathbf{f}_f \\ \mathbf{f}_c
    \end{bmatrix}.

Solving for :math:`\mathbf{u}_f` reduces to

.. math::

    \mathbf{K}_{ff}\, \mathbf{u}_f =
    \mathbf{f}_f - \mathbf{K}_{fc}\, \mathbf{u}_c,

and the constrained-DOF reactions follow as
:math:`\mathbf{r}_c = \mathbf{K}_{cf}\, \mathbf{u}_f +
\mathbf{K}_{cc}\, \mathbf{u}_c - \mathbf{f}_c`.

References
----------

* Bathe, K.-J. (2014) *Finite Element Procedures*, 2nd ed.,
  Watertown, MA, §4 (linear FE for solid mechanics).
* Zienkiewicz, O. C. and Taylor, R. L. (2013) *The Finite
  Element Method: Its Basis and Fundamentals*, 7th ed.,
  Butterworth-Heinemann, §2 (variational principles), §8
  (general elastic body equations).
* Hughes, T. J. R. (2000) *The Finite Element Method —
  Linear Static and Dynamic Finite Element Analysis*, Dover,
  §1 (linear elastostatics) + §3 (Galerkin discretisation).
* Cook, R. D., Malkus, D. S., Plesha, M. E., Witt, R. J. (2002)
  *Concepts and Applications of Finite Element Analysis*, 4th
  ed., Wiley, §4 (variational methods).
* Reddy, J. N. (2004) *An Introduction to Nonlinear Finite
  Element Analysis*, Oxford UP, §2 (weak-form derivation
  including dynamic and damping terms).