Define user problems

Fedoo can be extended at two complementary levels:

a weak form defines the equation assembled on the mesh;
a constitutive law defines the material response used by a mechanical weak form.

For a linear equation, a DiffOp expression and fedoo.WeakForm are generally sufficient. For a nonlinear or history-dependent equation, derive from WeakFormBase. For a mechanical material, MechanicalUMAT is the recommended extension point; derive directly from Mechanical3D only when full control over the Fedoo constitutive lifecycle is required.

Define a weak form from differential operators

A modeling space owns the unknown fields and creates their differential operators. The virtual property selects the test-function side of a weak term. Ordinary Python arithmetic combines these objects into a DiffOp expression.

The following example defines the weak form of the Poisson equation

\[\int_\Omega \nabla \delta u \mathbin{\cdot} \nabla u\,d\Omega = 0.\]

import fedoo as fd

space = fd.ModelingSpace("2D")
space.new_variable("U")

dU_dX = space.derivative("U", "X")
dU_dY = space.derivative("U", "Y")

equation = dU_dX.virtual * dU_dX + dU_dY.virtual * dU_dY
weakform = fd.WeakForm(equation)

mesh = fd.mesh.rectangle_mesh()
assembly = fd.Assembly.create(weakform, mesh)
problem = fd.problem.Linear(assembly)

problem.bc.add("Dirichlet", "left", "U", 0)
problem.bc.add("Dirichlet", "right", "U", 1)
problem.solve()

The same construction applies to systems of equations. Declare every scalar unknown with fedoo.ModelingSpace.new_variable(), group components with fedoo.ModelingSpace.new_vector() when useful, and add the weak terms. Mixed formulations may assign different element interpolations to different variables.

Complete examples are available in poison_eq.py and fluid_mechanics.py.

State-dependent weak forms

When coefficients or residual terms depend on the current solution, derive a class from WeakFormBase. Its principal methods are:

initialize(assembly, problem): Allocate values required by the formulation in assembly.sv.
update(assembly, problem): Evaluate the current fields and update the residual state before assembly.
get_weak_equation(assembly, problem): Return the linearized DiffOp, including the tangent and current residual terms.
set_start(assembly, problem): Commit the state at the beginning of a new time increment.
to_start(assembly, problem): Restore the beginning-of-increment state after a rejected increment.
reset(): Clear the complete problem history when necessary.

The problem may have problem.dtime == 0 during initialization. A custom weak form must not divide by dtime in that state. Time-dependent terms can be omitted until a positive increment is supplied. See navier_stokes.py for a complete nonlinear example.

Define a mechanical constitutive law

Use MechanicalUMAT for a conventional material-point update. It handles the Fedoo lifecycle, state allocation, temperature lookup, material local frames, finite-rotation increments, and the transformations of strains, stresses, and tangents. The callback only integrates the constitutive equation in the supplied corotational material coordinates.

Recommended: `MechanicalUMAT`

The callback has the following signature:

def umat(
    strain,
    dstrain,
    F0,
    F1,
    stress,
    DR,
    props,
    statev_start,
    time,
    dtime,
    wm_start,
    temperature,
    *,
    ndi,
    tangent_mode,
):
    ...
    return stress_new, statev_new, wm_new, tangent

Its arguments follow these conventions:

Argument	Meaning
`strain`	Total strain at the beginning of the increment, shape `(6, n_points)`.
`dstrain`	Current strain increment, shape `(6, n_points)`.
`F0`, `F1`	Beginning and current deformation gradients, shape `(3, 3, n_points)` in a finite-strain analysis. They are empty in a small-strain analysis.
`stress`	Stress at the beginning of the increment, shape `(6, n_points)`.
`DR`	Objective incremental rotation, shape `(3, 3, n_points)`.
`props`	Material properties stored as a Fortran-contiguous column.
`statev_start`	Private state at the beginning of the increment, shape `(n_statev, n_points)`.
`time`, `dtime`	Current problem time and time increment.
`wm_start`	Beginning-of-increment work variables, shape `(4, n_points)`.
`temperature`	Gauss-point temperatures, or `None` when no temperature field is available.
`ndi`	Number of direct components: 3 normally and 2 for the plane-stress interface.
`tangent_mode`	User-defined tangent selector forwarded unchanged by Fedoo.

Stress and strain use engineering Voigt notation. The component ordering, including the axisymmetric convention, is documented by Mechanical3D.

The callback returns the updated stress, state variables, work variables, and the material tangent. A tangent can have shape (6, 6) when it is uniform, or (6, 6, n_points) when it varies between material points. Do not modify statev_start or wm_start in place: they represent the committed state and must remain available if the global increment is rejected.

Initialization call

Fedoo calls the UMAT once during initialization with:

time == 0 and dtime == 0;
zero strain, strain increment, stress, and work arrays;
the assembly’s initial state-variable values, or zeros when none were supplied;
identity DR;
identity F0 and F1 for finite strain, or empty deformation-gradient arrays for small strain.

The dtime == 0 call is not a physical material increment. It must return the initial elastic tangent without advancing irreversible state, accumulating work, or evaluating expressions that divide by dtime. This tangent is stored as ElasticMatrix when use_elastic_tangent=True and is restored at the beginning of each increment for the elastic predictor.

Initial state variables

State-variable initial conditions belong to an assembly rather than to the material, allowing one material instance to be reused by several problems with different histories. Set labeled values after creating the assembly and before initializing the problem:

import fedoo as fd

material = MechanicalUMAT(
    umat,
    props,
    n_statev=7,
    statev_label={"P": 0, "EP": slice(1, 7)},
)
weakform = fd.weakform.StressEquilibrium(material)
assembly = fd.Assembly.create(weakform, mesh)

material.set_initial_statev(assembly, "P", 0.02)
material.set_initial_statev(assembly, "EP", initial_plastic_strain)

A scalar is broadcast over the selected components and all Gauss points. A component vector is broadcast over the Gauss points, while a complete field has shape (n_components, assembly.n_gauss_points). For a scalar label, a one-dimensional array of length assembly.n_gauss_points supplies spatially varying values. Node and element fields are not converted implicitly.

Advanced users may instead assign the complete assembly.sv["Statev"] array manually before initialization. Its shape must be exactly (material.n_statev, assembly.n_gauss_points). In either case, the initial state is passed to the dtime == 0 UMAT call, so it may affect the initial tangent.

Calling assembly.reset() clears runtime state, including these initial values. Call set_initial_statev again before reinitializing a reset assembly.

Minimal elastic callback

This small 3D example illustrates the callback contract. Production laws should also validate their parameters and supported dimensional assumptions.

import numpy as np
from fedoo.core import MechanicalUMAT


def elastic_umat(
    strain, dstrain, F0, F1, stress, DR, props,
    statev_start, time, dtime, wm_start, temperature,
    *, ndi, tangent_mode,
):
    E, nu = props[:, 0]
    shear = E / (2 * (1 + nu))
    lame = E * nu / ((1 + nu) * (1 - 2 * nu))

    tangent = np.zeros((6, 6))
    tangent[:3, :3] = lame
    tangent[0, 0] += 2 * shear
    tangent[1, 1] += 2 * shear
    tangent[2, 2] += 2 * shear
    tangent[3, 3] = shear
    tangent[4, 4] = shear
    tangent[5, 5] = shear

    stress_new = stress + tangent @ dstrain
    return (
        stress_new,
        statev_start.copy(),
        wm_start.copy(),
        tangent,
    )


material = MechanicalUMAT(
    elastic_umat,
    props=[210_000.0, 0.3],
    props_label={"E": 0, "nu": 1},
    is_isotropic=True,
)

The constructor also accepts:

n_statev and statev_label: Allocate private history and expose selected components to post-processing.
density: Define the mass density for dynamic mechanical problems.
tangent_from_F: Mark a finite-strain law whose tangent is derived from deformation-gradient kinematics instead of the default corotational strain measure. Such a law requires a geometrically nonlinear formulation; F0 and F1 remain empty in small strain.
use_elastic_tangent: Control whether the initialization tangent is restored at the beginning of each increment.
tangent_mode: Select the tangent requested from the callback.
is_isotropic: Avoid unnecessary frame transformations for an invariant response.

Local frames and finite rotations

Assigning a frame with material.set_local_frame(...) defines the initial material coordinates. MechanicalUMAT transforms the strain, deformation gradient, starting stress, and incremental rotation into those coordinates, then transforms the returned stress and tangent back to the global frame. Uniform, nodal, elemental, and Gauss-point material frames are supported as described in Constitutive Law (fedoo.constitutivelaw).

Private state variables are opaque to Fedoo because their tensorial meaning is law-dependent. The callback must apply DR to any tensor-valued history that requires objective transport. Scalar history variables require no rotation.

The objective-rate formulation is selected on fedoo.weakform.StressEquilibrium through weakform.corate. A UMAT subclass with a strict requirement can set, for example:

class MyUMAT(MechanicalUMAT):
    required_corate = ("log_r", "log_r_inc")

Fedoo then validates the selected formulation when the material is initialized.

Advanced: derive directly from `Mechanical3D`

Deriving from Mechanical3D provides complete control but also transfers the constitutive lifecycle to the implementation. This route is appropriate when the response cannot be expressed by the UMAT contract, needs custom assembly state, or uses a specialized tangent representation.

A derived law is generally responsible for:

calling super().__init__(name=name, density=density);
declaring is_isotropic correctly;
allocating its state and initial assembly.sv["TangentMatrix"] in initialize;
reading assembly.sv["Strain"] and assembly.sv["DStrain"] and writing assembly.sv["Stress"] and assembly.sv["TangentMatrix"] in update;
implementing commit, rollback, and reset behavior through set_start, to_start, and reset when it owns history variables;
handling plane stress explicitly when that assumption is supported;
using the frame-transformation helpers supplied by Mechanical3D when the law is anisotropic;
setting _Lt_from_F only when the tangent follows deformation-gradient rather than strain-increment kinematics.

assembly.sv_start is a shallow beginning-of-increment snapshot. Do not modify arrays shared with it in place during a trial iteration. Replace trial arrays in assembly.sv or make explicit copies so that to_start can restore the committed state.

Direct Mechanical3D laws must also tolerate initialization before a positive time increment exists. They should always provide the initial tangent required to assemble the first global system. The implementations of fedoo.constitutivelaw.ElasticAnisotropic and the pedagogical fedoo.constitutivelaw.ElastoPlasticity provide useful references; the latter uses MechanicalUMAT for its lifecycle.