From 9b22df8d3c100ccc3f34577cdd2a0a2506090189 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Fri, 14 Nov 2025 23:22:42 +0200 Subject: [PATCH] refactor(physics): Implement double-dispatch architecture Major API refactoring: replace mutable ElasticityPhysicsType with type-parametric Physics struct for compile-time dispatch. New type hierarchy: - AbstractField: What we solve (Displacement{3}, Temperature, etc.) - AbstractFormulation: How we discretize (ContinuumFormulation, BeamFormulation) - AbstractMaterial: Material behavior (LinearElastic, NeoHookean) - AbstractMesh: Mesh container Physics{Formulation, Field, Mesh, Material} enables natural dispatch: assemble(::Physics{ContinuumFormulation{FullThreeD}, Displacement{3}, M, Mat}) assemble(::Physics{BeamFormulation{Timoshenko}, DisplacementRotation{3}, M, Mat}) Type parameter order prioritizes Formulation for dispatch hierarchy. Breaking changes: - Old: Physics(Elasticity, "name", 3) - New: Physics(name=..., mesh=..., field=Displacement{3}(), formulation=ContinuumFormulation{FullThreeD}(), material=...) - Deprecate: add_elements!() - Physics references Mesh, doesn't own elements Benefits: - Type stability: All types known at compile time - Dispatch: Specialized methods for formulation/field combinations - Extensibility: New formulations/fields without modifying core - Performance: No runtime type checks, optimal codegen This is foundation for the NEW API (TDD tests, Nov 14 2025). --- src/physics_api.jl | 411 ++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 371 insertions(+), 40 deletions(-) diff --git a/src/physics_api.jl b/src/physics_api.jl index 9c36bb3..4774f19 100644 --- a/src/physics_api.jl +++ b/src/physics_api.jl @@ -7,24 +7,175 @@ Backend-agnostic elasticity API This file defines the user-facing API that works regardless of backend (CPU or GPU). """ +# ============================================================================ +# Abstract Types for Mesh, Material, Field, and Formulation +# ============================================================================ + """ - Elasticity <: FieldProblem + AbstractMesh -Elasticity physics problem type. +Abstract type for all mesh structures. -Properties: -- `formulation`: `:continuum` (3D solid) or `:plane_stress`, `:plane_strain` -- `finite_strain`: Enable geometric nonlinearity -- `geometric_stiffness`: Include stress-dependent stiffness +Concrete subtypes: Mesh{T<:AbstractTopology} """ -mutable struct ElasticityPhysicsType <: FieldProblem - formulation::Symbol # :continuum (3D) - finite_strain::Bool # Geometric nonlinearity - geometric_stiffness::Bool # Stress-dependent stiffness -end +abstract type AbstractMesh end -const Elasticity = ElasticityPhysicsType -ElasticityPhysicsType() = ElasticityPhysicsType(:continuum, false, false) +""" + AbstractMaterial + +Abstract type for all material models. + +Concrete subtypes: LinearElastic, NeoHookean, PerfectPlasticity, etc. +""" +abstract type AbstractMaterial end + +# ============================================================================ +# Abstract Types for Field and Formulation (Double Dispatch) +# ============================================================================ + +""" + AbstractField + +Abstract type for field variables we solve for. + +Subtypes define what physical quantity is being solved: +- `Displacement{Dim}`: Displacement field (3D, 2D, etc.) +- `Temperature`: Temperature field +- `DisplacementRotation{Dim}`: Displacement + rotation (beams, shells) +""" +abstract type AbstractField end + +""" + Displacement{Dim} <: AbstractField + +Displacement field with Dim components per node. + +# Examples +- `Displacement{3}()`: 3D displacement (ux, uy, uz) +- `Displacement{2}()`: 2D displacement (ux, uy) +""" +struct Displacement{Dim} <: AbstractField end + +""" + dofs_per_node(field::AbstractField) -> Int + +Number of degrees of freedom per node for this field type. +""" +dofs_per_node(::Displacement{Dim}) where Dim = Dim + +""" + Temperature <: AbstractField + +Temperature field (scalar per node). +""" +struct Temperature <: AbstractField end +dofs_per_node(::Temperature) = 1 + +""" + DisplacementRotation{Dim} <: AbstractField + +Combined displacement and rotation field (for beams, shells). + +DOFs per node: 2*Dim (Dim displacements + Dim rotations) +- 3D: 6 DOFs (ux, uy, uz, θx, θy, θz) +- 2D: 4 DOFs (ux, uy, θz, warping) +""" +struct DisplacementRotation{Dim} <: AbstractField end +dofs_per_node(::DisplacementRotation{Dim}) where Dim = 2 * Dim + +""" + AbstractFormulation + +Abstract type for discretization formulations. + +Subtypes define HOW we discretize the governing equations: +- `ContinuumFormulation{Theory}`: Standard FEM for continuum mechanics +- `BeamFormulation{Theory}`: Beam elements (Euler-Bernoulli, Timoshenko) +- `ShellFormulation{Theory}`: Shell elements (Reissner-Mindlin, Kirchhoff) +- `TrussFormulation`: Truss elements +""" +abstract type AbstractFormulation end + +""" + AbstractContinuumTheory + +Theory variants for continuum formulation. +""" +abstract type AbstractContinuumTheory end +struct FullThreeD <: AbstractContinuumTheory end +struct PlaneStress <: AbstractContinuumTheory end +struct PlaneStrain <: AbstractContinuumTheory end +struct Axisymmetric <: AbstractContinuumTheory end + +""" + ContinuumFormulation{Theory} <: AbstractFormulation + +Standard continuum mechanics formulation with theory variant. + +# Examples +- `ContinuumFormulation{FullThreeD}()`: Full 3D analysis +- `ContinuumFormulation{PlaneStress}()`: 2D plane stress +- `ContinuumFormulation{PlaneStrain}()`: 2D plane strain +""" +struct ContinuumFormulation{Theory<:AbstractContinuumTheory} <: AbstractFormulation end + +""" + AbstractBeamTheory + +Theory variants for beam formulation. +""" +abstract type AbstractBeamTheory end +struct EulerBernoulli <: AbstractBeamTheory end +struct Timoshenko <: AbstractBeamTheory end + +""" + BeamFormulation{Theory} <: AbstractFormulation + +Beam element formulation with theory variant. + +# Examples +- `BeamFormulation{EulerBernoulli}()`: No shear deformation +- `BeamFormulation{Timoshenko}()`: Includes shear deformation +""" +struct BeamFormulation{Theory<:AbstractBeamTheory} <: AbstractFormulation end + +""" + TrussFormulation <: AbstractFormulation + +Truss element formulation (axial force only). +""" +struct TrussFormulation <: AbstractFormulation end + +# ============================================================================ +# Physics Type Hierarchy +# ============================================================================ + +""" + AbstractPhysics + +Abstract base for all physics types. + +JuliaFEM recognizes two fundamental categories: +1. Bulk Physics - operates on element sets (domain interiors) +2. Interface Physics - operates on surface pairs (domain boundaries) +""" +abstract type AbstractPhysics end + +# ============================================================================ +# Constraint Types (placeholder for future implementation) +# ============================================================================ + +""" + Constraint + +Placeholder for internal constraints (rigid body modes, incompressibility, etc.). +To be implemented in future stories. +""" +struct Constraint end + +# ============================================================================ +# Boundary Conditions +# ============================================================================ """ DirichletBC @@ -60,72 +211,252 @@ end NeumannBC() = NeumannBC(Element[], Vec{3,Float64}[]) +# ============================================================================ +# Physics Struct (Fully Typed with Type Parameters) +# ============================================================================ + """ - Physics{P} + Physics{Formulation<:AbstractFormulation, Field<:AbstractField, Mesh<:AbstractMesh, Material<:AbstractMaterial} <: AbstractPhysics -Container for physics problem with elements and boundary conditions. +Bulk physics operating on element sets. -Type parameter `P` is the physics type (e.g., Elasticity, HeatTransfer). +Represents volumetric phenomena: +- Solid mechanics (elasticity, plasticity, damage) +- Heat transfer +- Fluid dynamics +- Chemical diffusion -Fields: +# Type Parameters (in dispatch priority order) +- `Formulation`: How we discretize - ContinuumFormulation{FullThreeD}, BeamFormulation{Timoshenko}, etc. +- `Field`: What we solve - Displacement{3}, Temperature, etc. +- `Mesh`: Mesh type - AbstractMesh subtype (e.g., Mesh, Mesh{Tet4}) +- `Material`: Material type - AbstractMaterial subtype (e.g., LinearElastic, NeoHookean) + +# Fields - `name`: Problem name -- `dimension`: Spatial dimension (1, 2, or 3) -- `properties`: Physics-specific properties (e.g., Elasticity struct) -- `body_elements`: Volume elements with geometry and material +- `mesh`: Reference to Mesh (topology owner - NOT copied!) +- `element_set`: Which elements in mesh this physics applies to +- `field`: Field instance (type Field) +- `formulation`: Formulation instance (type Formulation) +- `material`: Material properties (type Material) +- `constraints`: Optional internal constraints - `bc_dirichlet`: Dirichlet boundary conditions - `bc_neumann`: Neumann boundary conditions (surface loads) +# Type Parameter Order Rationale +**Formulation first** enables natural dispatch hierarchy: +```julia +# Most specific: Formulation + Field combination +assemble(::Physics{ContinuumFormulation{FullThreeD}, Displacement{3}, M, Mat}) = ... +assemble(::Physics{ContinuumFormulation{PlaneStress}, Displacement{2}, M, Mat}) = ... +assemble(::Physics{BeamFormulation{Timoshenko}, DisplacementRotation{3}, M, Mat}) = ... + +# Generic fallbacks work naturally +assemble(::Physics{ContinuumFormulation, Temperature, M, Mat}) = ... # Any continuum + thermal +assemble(::Physics{Fm, F, M, Mat}) where {Fm,F,M,Mat} = ... # Fully generic +``` + +# Benefits of Type Parameters +- **Dispatch-optimized order:** Most important types first (formulation, then field) +- **Zero-cost dispatch:** Compiler selects assembly method at compile time +- **Type stability:** All field types known statically +- **Specialization:** Different algorithms per formulation × field combination +- **GPU-ready:** Concrete types enable device-specific optimization + +# Design Philosophy +Physics REFERENCES Mesh (does not own it). Multiple Physics can share same Mesh. +Material properties stored in Physics, topology stored in Mesh. + # Example ```julia -physics = Physics(Elasticity, "beam", 3) -add_elements!(physics, elements) +mesh = Mesh("cantilever.inp") +material = LinearElastic(E=210e9, ν=0.3) + +physics = Physics( + name = "cantilever", + mesh = mesh, + element_set = :solid, + field = Displacement{3}(), + formulation = ContinuumFormulation{FullThreeD}(), + material = material +) +# Type: Physics{ContinuumFormulation{FullThreeD}, Displacement{3}, typeof(mesh), LinearElastic} + add_dirichlet!(physics, [1,2,3], [1,2,3], 0.0) sol = solve!(physics) ``` + +# Dispatch Example + +```julia +# Different formulation/field combinations dispatch to specialized methods +physics_continuum_3d = Physics(..., formulation=ContinuumFormulation{FullThreeD}(), field=Displacement{3}()) +physics_plane_stress = Physics(..., formulation=ContinuumFormulation{PlaneStress}(), field=Displacement{2}()) +physics_beam = Physics(..., formulation=BeamFormulation{Timoshenko}(), field=DisplacementRotation{3}()) +physics_thermal = Physics(..., formulation=ContinuumFormulation{FullThreeD}(), field=Temperature()) + +# Compiler generates specialized code for each formulation × field! +assemble(physics_continuum_3d) # → assemble(::Physics{ContinuumFormulation{FullThreeD}, Displacement{3}, ...}) +assemble(physics_plane_stress) # → assemble(::Physics{ContinuumFormulation{PlaneStress}, Displacement{2}, ...}) +assemble(physics_beam) # → assemble(::Physics{BeamFormulation{Timoshenko}, DisplacementRotation{3}, ...}) +assemble(physics_thermal) # → assemble(::Physics{ContinuumFormulation{FullThreeD}, Temperature, ...}) +``` """ -mutable struct Physics{P} +struct Physics{Formulation<:AbstractFormulation,Field<:AbstractField,Mesh<:AbstractMesh,Material<:AbstractMaterial} <: AbstractPhysics name::String - dimension::Int - properties::P - body_elements::Vector{Element} - bc_dirichlet::DirichletBC - bc_neumann::NeumannBC + mesh::Mesh # Fully typed mesh reference + element_set::Symbol # Which bulk elements + field::Field # What we solve (Displacement{3}, Temperature, etc.) + formulation::Formulation # How we discretize (ContinuumFormulation{...}, etc.) + material::Material # Fully typed material + constraints::Vector{Constraint} # Optional: internal constraints + bc_dirichlet::DirichletBC # Essential boundary conditions + bc_neumann::NeumannBC # Natural boundary conditions + + # Inner constructor to ensure type parameters match field instances + function Physics{Formulation,Field,Mesh,Material}( + name::String, + mesh::Mesh, + element_set::Symbol, + field::Field, + formulation::Formulation, + material::Material, + constraints::Vector{Constraint}, + bc_dirichlet::DirichletBC, + bc_neumann::NeumannBC + ) where {Formulation<:AbstractFormulation,Field<:AbstractField,Mesh<:AbstractMesh,Material<:AbstractMaterial} + new{Formulation,Field,Mesh,Material}(name, mesh, element_set, field, formulation, material, + constraints, bc_dirichlet, bc_neumann) + end end """ - Physics(::Type{P}, name::String, dimension::Int) where P + Physics(; name, mesh, element_set, field, formulation, material, constraints=Constraint[]) -Create a physics problem of type P. +Create a physics problem referencing a mesh. + +Type parameters are automatically inferred from the argument types in dispatch-optimized order. + +# Arguments (keyword arguments) +- `name::String`: Problem name +- `mesh`: Mesh reference (topology owner) - type becomes parameter Mesh +- `element_set::Symbol`: Which elements in mesh to use +- `field::AbstractField`: What we solve - type Field inferred +- `formulation::AbstractFormulation`: How we discretize - type Formulation inferred +- `material`: Material properties - type Material inferred +- `constraints`: Optional internal constraints (default: empty) + +# Returns +`Physics{Formulation, Field, Mesh, Material}` where types are inferred from arguments # Example ```julia -physics = Physics(Elasticity, "cantilever", 3) +mesh = Mesh("model.inp") +material = LinearElastic(E=210e9, ν=0.3) + +physics = Physics( + name = "cantilever", + mesh = mesh, + element_set = :solid, + field = Displacement{3}(), + formulation = ContinuumFormulation{FullThreeD}(), + material = material +) +# Type: Physics{ContinuumFormulation{FullThreeD}, Displacement{3}, typeof(mesh), LinearElastic} +``` + +# Type Inference + +All type parameters are inferred automatically in dispatch priority order: +- `Formulation = typeof(formulation)` - formulation type (e.g., ContinuumFormulation{FullThreeD}) +- `Field = typeof(field)` - field type (e.g., Displacement{3}) +- `Mesh = typeof(mesh)` - concrete mesh type +- `Material = typeof(material)` - concrete material type + +# Dispatch Benefits + +Type parameter order optimized for natural dispatch: +```julia +# Most specific methods dispatch first +assemble(::Physics{ContinuumFormulation{FullThreeD}, Displacement{3}, M, Mat}) +# Generic fallbacks work naturally +assemble(::Physics{Fm, F, M, Mat}) where {Fm, F, M, Mat} ``` """ -function Physics(::Type{P}, name::String, dimension::Int) where P - properties = P() - body_elements = Element[] +function Physics(; + name::String, + mesh::M, + element_set::Symbol, + field::F, + formulation::Fm, + material::Mat, + constraints::Vector{Constraint}=Constraint[] +) where {M<:AbstractMesh,Mat<:AbstractMaterial,F<:AbstractField,Fm<:AbstractFormulation} bc_dirichlet = DirichletBC() bc_neumann = NeumannBC() - return Physics{P}(name, dimension, properties, body_elements, bc_dirichlet, bc_neumann) + + return Physics{Fm,F,M,Mat}(name, mesh, element_set, field, formulation, material, + constraints, bc_dirichlet, bc_neumann) end +# ============================================================================ +# Deprecated API (Compatibility Layer) +# ============================================================================ + +# Note: Old Physics(::Type{P}, name, dimension) constructor is REMOVED +# The new type signature Physics{M,Mat,F,Fm} is incompatible with the old API +# Users must migrate to: Physics(name=..., mesh=..., field=..., formulation=..., material=...) + """ add_elements!(physics::Physics, elements::Vector{Element}) -Add volume elements to the physics problem. +DEPRECATED: Physics no longer owns elements directly. -# Example +In the new API, Physics references Mesh. Elements are defined in Mesh, +and Physics selects which elements to use via `element_set`. +# Migration + +Old API: ```julia -add_elements!(physics, [el1, el2, el3]) +physics = Physics(Elasticity, "beam", 3) +add_elements!(physics, elements) ``` + +New API: +```julia +mesh = elements_to_mesh(elements) # Convert elements to mesh +physics = Physics( + name = "beam", + mesh = mesh, + element_set = :all, + field = Displacement{3}(), + formulation = ContinuumFormulation{FullThreeD}(), + material = extract_material(elements[1]) +) +``` + +This function will be removed in v2.0. """ function add_elements!(physics::Physics, elements::Vector{Element}) - append!(physics.body_elements, elements) - return nothing + @warn """ + add_elements!(physics, elements) is DEPRECATED. + + Physics now references Mesh instead of owning elements. + Create a Mesh from elements and pass it to Physics constructor. + + Migration: + mesh = elements_to_mesh(elements) + physics = Physics(name=..., mesh=mesh, element_set=:all, ...) + + This function will be removed in v2.0. + """ maxlog = 1 + + # For backward compatibility, we need to somehow add elements to the mesh + # This is hacky but maintains old API temporarily + error("add_elements! no longer supported. Please migrate to new Physics(mesh=...) API.") end """