diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..e86d952 --- /dev/null +++ b/.gitignore @@ -0,0 +1,5 @@ +# ignore cached bytecode-compiled python scene files +*__pycache__/ + +# a run's own output +verification/results/ diff --git a/README.md b/README.md index 25d0262..497734e 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,45 @@ # Verification-Validation Initially funded by TIRREX, this repository proposes a framework for verification and validation. It relies on both manufactured solutions and comparison with alternative simulation technologies. This work therefore proves the reliability and compliance of the SOFA models & algorithms. + +## Verification + +**Verification** uses the method of manufactured solutions: an exact displacement is chosen, the +loads that produce it are derived, and SOFA solves the problem on a sequence of refined meshes. +The error must shrink at the rate the discretisation predicts. Covered force fields: +- `LinearSmallStrainFEMForceField` + +To get started, read [`verification/docs/IntroductionMMS.ipynb`](verification/docs/IntroductionMMS.ipynb): it presents the method and builds one verification study step by step, from the manufactured solution to the input file the suite runs. + +The following force field components are currently verified within this repository: +- `LinearSmallStrainFEMForceField` [`verification/docs/LinearSmallStrainFEMForceField.ipynb`](verification/docs/LinearSmallStrainFEMForceField.ipynb) + +## Layout + +``` +common/ shared code: geometry and SOFA helpers +verification/ the verification suite + cases//**.json one test case per file +``` + +## Requirements + +- SOFA with SofaPython3 +- numpy, sympy, matplotlib + +## Running + +```bash +python verification/run.py # one test case +python verification/run.py --all # every test case, then one summary line per case +python verification/run.py --clean # delete previous results +``` + +Results are written to `verification/results/`. + +To inspect one mesh of a case in the SOFA GUI, the runner also acts as a SOFA scene file: + +```bash +runSofa -l SofaPython3 verification/run.py --argv --argv 0 +``` + +The first argument is the path to the .json case file and the second is the mesh refinement level to be built. diff --git a/common/__init__.py b/common/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/common/geometry/__init__.py b/common/geometry/__init__.py new file mode 100644 index 0000000..f82eb80 --- /dev/null +++ b/common/geometry/__init__.py @@ -0,0 +1,3 @@ +from .geometry import Geometry +from .grid import Grid, Bar1D, Beam2D, Beam3D +from .tire import Tire, Tire2D, Tire3D diff --git a/common/geometry/geometry.py b/common/geometry/geometry.py new file mode 100644 index 0000000..46e44ef --- /dev/null +++ b/common/geometry/geometry.py @@ -0,0 +1,116 @@ +"""Tool-agnostic geometry descriptions for the VnV suite.""" + +from abc import ABC, abstractmethod + + +def _region_indices_from_mesh(mesh_path): + """Node indices tagged by each physical group of a Gmsh .msh file, keyed by group name.""" + import meshio + + mesh = meshio.read(mesh_path) + tag_names = {tag: name for name, (tag, _dim) in mesh.field_data.items()} + groups = {} + for block, tags in zip(mesh.cells, mesh.cell_data.get("gmsh:physical", [])): + for element, tag in zip(block.data, tags): + name = tag_names.get(tag) + if name is not None: + groups.setdefault(name, set()).update(element.tolist()) + return {name: sorted(indices) for name, indices in groups.items()} + + +class Geometry(ABC): + """A shape parameterization plus its named boundary regions.""" + + dim: int # topological dimension of the shape (and of the elements meshing it) + PARAMETER_ORDER: tuple # the family's parameters, ordered by the dimension that introduces each + + def __init__(self, *, dim=None, spatialDimensions=None, meshPath=None, **parameters): + """Build this shape, taking each name in its `PARAMETER_NAMES` as a keyword. + + e.g. ``Beam3D(length=2.0, width=1.0, height=1.0)``. Pass `dim` only for a family base, + which names no dimension of its own: ``Grid(dim=3, length=2.0, width=1.0, height=1.0)`` + builds that same box. + """ + stated = getattr(type(self), "dim", None) + if stated is None and dim is None: + raise TypeError(f"{type(self).__name__} states no dim: pass dim= to build one, or " + f"name the geometry that states it") + if stated is not None and dim is not None and dim != stated: + raise TypeError(f"{type(self).__name__} is {stated}D and cannot be built with dim={dim}") + self.dim = stated if dim is None else dim + + if set(parameters) != set(self.PARAMETER_NAMES): + raise TypeError(f"{type(self).__name__} takes {list(self.PARAMETER_NAMES)}, " + f"got {sorted(parameters)}") + + for name in self.PARAMETER_NAMES: + setattr(self, name, parameters[name]) + self._set_spatial_dimensions(spatialDimensions) + + if meshPath is not None: + self._regions = {name: (lambda nodes, indices=indices: indices) + for name, indices in _region_indices_from_mesh(meshPath).items()} + else: + self._regions = {name: (lambda nodes, predicate=predicate: + [i for i, p in enumerate(nodes) if predicate(p)]) + for name, predicate in self._analytic_predicates.items()} + + @property + def PARAMETER_NAMES(self) -> tuple: + """The parameters this shape takes: the first `dim` of its family's. + + Each dimension adds one, so a member takes as many as it has dimensions: of Grid's + ("length", "width", "height"), Bar1D takes ("length",) and Beam2D ("length", "width"). + A family that does not grow this way sets `PARAMETER_NAMES` itself instead. + """ + return self.PARAMETER_ORDER[:self.dim] + + @property + def parameters(self) -> list: + """This shape's defining parameters, read off the attributes named in `PARAMETER_NAMES`. + + e.g. Beam2D(length=2.0, width=0.5).parameters == [2.0, 0.5] + """ + return [getattr(self, name) for name in self.PARAMETER_NAMES] + + @property + def named_parameters(self) -> dict: + """Parameters keyed by name: how a manufactured solution asks for one of them. + + e.g. Beam2D(length=2.0, width=0.5).named_parameters == {"length": 2.0, "width": 0.5} + """ + return dict(zip(self.PARAMETER_NAMES, self.parameters)) + + def _set_spatial_dimensions(self, spatial_dimensions): + """Set how many dimensions of space this shape is embedded in; defaults to its own.""" + if spatial_dimensions is None: + spatial_dimensions = self.dim + + if spatial_dimensions < self.dim: + raise ValueError( + f"{type(self).__name__} is {self.dim}D and cannot be embedded in a " + f"{spatial_dimensions}D space: the embedding space must be at least as " + f"large as the shape." + ) + if spatial_dimensions > 3: + raise ValueError( + f"{type(self).__name__}: spatial_dimensions cannot exceed 3, " + f"got {spatial_dimensions}." + ) + + self.spatial_dimensions = spatial_dimensions + + @property + def region_names(self): + """This geometry's named boundary regions.""" + return self._regions.keys() + + def region_indices(self, region, nodes) -> list: + """Indices into `nodes` belonging to `region`.""" + return self._regions[region](nodes) + + @property + @abstractmethod + def _analytic_predicates(self) -> dict: + """Named boundary regions as ``{name: predicate(point) -> bool}``, from this geometry's + own shape parameters.""" diff --git a/common/geometry/grid.py b/common/geometry/grid.py new file mode 100644 index 0000000..d441498 --- /dev/null +++ b/common/geometry/grid.py @@ -0,0 +1,36 @@ +"""Rectilinear structured grid geometry in 1 to 3 dimensions: bars and beams.""" + +import numpy as np + +from .geometry import Geometry + + +class Grid(Geometry): + """A rectilinear structured grid, axis i spanning 0..parameters[i]. Always starts at the origin.""" + + PARAMETER_ORDER = ("length", "width", "height") + AXIS_NAMES = (("left", "right"), ("bottom", "top"), ("front", "back")) + + @property + def _analytic_predicates(self) -> dict: + predicates = {} + for axis, (low, high) in enumerate(self.AXIS_NAMES[:self.dim]): + extent = self.parameters[axis] + predicates[low] = lambda p, axis=axis: np.isclose(p[axis], 0.0) + predicates[high] = lambda p, axis=axis, extent=extent: np.isclose(p[axis], extent) + return predicates + + +class Bar1D(Grid): + """A 1D bar occupying ``[0, length]`` along x.""" + dim = 1 + + +class Beam2D(Grid): + """A 2D rectangle occupying ``[0, length] x [0, width]`` in the xy-plane.""" + dim = 2 + + +class Beam3D(Grid): + """A 3D box occupying ``[0, length] x [0, width] x [0, height]``.""" + dim = 3 diff --git a/common/geometry/tire.py b/common/geometry/tire.py new file mode 100644 index 0000000..2eef1df --- /dev/null +++ b/common/geometry/tire.py @@ -0,0 +1,32 @@ +"""Tire geometry (Chamberland et al., 2010, Fig. 1).""" + +import numpy as np + +from .geometry import Geometry + + +class Tire(Geometry): + """An annulus in the xy-plane, centered on the origin, extruded along z when it is 3D.""" + + PARAMETER_ORDER = ("inner_radius", "outer_radius", "extrusion_length") + + @property + def _analytic_predicates(self) -> dict: + predicates = { + "inner": lambda p: np.isclose(np.hypot(p[0], p[1]), self.inner_radius), + "outer": lambda p: np.isclose(np.hypot(p[0], p[1]), self.outer_radius), + } + if self.dim == 3: + predicates["front"] = lambda p: np.isclose(p[2], 0.0) + predicates["back"] = lambda p: np.isclose(p[2], self.extrusion_length) + return predicates + + +class Tire2D(Tire): + """An annulus occupying ``inner_radius <= r <= outer_radius`` of the xy-plane.""" + dim = 2 + + +class Tire3D(Tire): + """A Tire2D extruded along z from 0 to extrusion_length.""" + dim = 3 diff --git a/common/sofa/__init__.py b/common/sofa/__init__.py new file mode 100644 index 0000000..3175991 --- /dev/null +++ b/common/sofa/__init__.py @@ -0,0 +1 @@ +from .scene import SofaScene diff --git a/common/sofa/controllers.py b/common/sofa/controllers.py new file mode 100644 index 0000000..4bd21e3 --- /dev/null +++ b/common/sofa/controllers.py @@ -0,0 +1,47 @@ +"""Post-init SOFA controllers.""" + +import Sofa.Core + + +class RegionClamp(Sofa.Core.Controller): + """Builds one PartialFixedProjectiveConstraint per boundary region, then applies Dirichlet + boundary conditions once the simulation has initialized. + + For each region, the nodes belonging to it are found and handed to its constraint as the indices + it fixes. If a prescribed displacement is given, those nodes are also moved from rest to rest + + u in the constraint's fixed directions before the constraint clamps them there. rest_position + itself is never changed. + """ + + def __init__(self, geometry, dofs, node, prescribe_displacement_on, vec_type, + displacement=None, *args, **kwargs): + super().__init__(*args, **kwargs) + self.geometry = geometry + self.dofs = dofs + self.displacement = displacement + + # Fix dofs partially as determined by MMS `prescribe_displacement_on` + self.region_constraints = [] + for region, fixed_directions in prescribe_displacement_on.items(): + # Pad the array with zeroes for embedded dimensions: e.g. [1, 0] -> [1, 0, 0] + padded = tuple(fixed_directions) + (0,) * (geometry.spatial_dimensions - len(fixed_directions)) + constraint = node.addObject('PartialFixedProjectiveConstraint', name='clamp' + region, + template=vec_type, fixedDirections=list(padded)) + self.region_constraints.append((region, padded, constraint)) + + def onSimulationInitDoneEvent(self, event): + rest_positions = self.dofs.rest_position.array() + + # Apply non-zero Dirichlet BC: Move each fixed direction to rest + u; rest_position is untouched. + with self.dofs.position.writeableArray() as pos: + for region, fixed_directions, constraint in self.region_constraints: + indices = self.geometry.region_indices(region, rest_positions) + + if self.displacement is not None and indices: + displacement = self.displacement(rest_positions[indices]) + + for axis, fixed in enumerate(fixed_directions): + if fixed: + pos[indices, axis] = rest_positions[indices, axis] + displacement[:, axis] + + constraint.indices.value = indices diff --git a/common/sofa/conventions.py b/common/sofa/conventions.py new file mode 100644 index 0000000..d14acd8 --- /dev/null +++ b/common/sofa/conventions.py @@ -0,0 +1,38 @@ +"""Central map from tool-agnostic vocabulary to SOFA components.""" + +from dataclasses import dataclass + +# Spatial dimension (DOF/embedding space) -> MechanicalObject / force-field vector template. +VEC_DIM = {1: "Vec1d", 2: "Vec2d", 3: "Vec3d"} + +@dataclass(frozen=True) +class ElementKind: + """Everything the scene-building layer needs to know about one type of finite element.""" + + dim: int # topological dimension of this element + container: str # topology container class + data_name: str # Data field this kind's own container addresses its elements by + cpp: str # SOFA geometry name, for compound templates ("Vec3d,Hexahedron") + facet_kind: str | None # the kind of this element's own facets, None if it has none + grid_mapping: str | None # topological mapping generating it from the grid, None if grid-native + boundary_mappings: tuple # (element kind, mapping) steps from this element to its boundary + +ELEMENTS = { + "edge": ElementKind(dim=1, container="EdgeSetTopologyContainer", data_name="edges", + cpp="Edge", facet_kind=None, grid_mapping=None, + boundary_mappings=()), + "tri": ElementKind(dim=2, container="TriangleSetTopologyContainer", data_name="triangles", + cpp="Triangle", facet_kind="edge", grid_mapping="Quad2TriangleTopologicalMapping", + boundary_mappings=(("edge", "Triangle2EdgeTopologicalMapping"),)), + "quad": ElementKind(dim=2, container="QuadSetTopologyContainer", data_name="quads", + cpp="Quad", facet_kind="edge", grid_mapping=None, + # To arrive at a Quad2Edge map -> triangulate quads and then Tri2Edge + boundary_mappings=(("tri", "Quad2TriangleTopologicalMapping"), + ("edge", "Triangle2EdgeTopologicalMapping"))), + "tet": ElementKind(dim=3, container="TetrahedronSetTopologyContainer", data_name="tetrahedra", + cpp="Tetrahedron", facet_kind="tri", grid_mapping="Hexa2TetraTopologicalMapping", + boundary_mappings=(("tri", "Tetra2TriangleTopologicalMapping"),)), + "hexa": ElementKind(dim=3, container="HexahedronSetTopologyContainer", data_name="hexahedra", + cpp="Hexahedron", facet_kind="quad", grid_mapping=None, + boundary_mappings=(("quad", "Hexa2QuadTopologicalMapping"),)), +} diff --git a/common/sofa/prefabs/__init__.py b/common/sofa/prefabs/__init__.py new file mode 100644 index 0000000..8ea7478 --- /dev/null +++ b/common/sofa/prefabs/__init__.py @@ -0,0 +1,8 @@ +from .base import ScenePrefab +from .grid import ElasticGrid, ElasticBar, ElasticBeam +from .tire import ElasticTire +from ...geometry import Grid, Bar1D, Beam2D, Beam3D, Tire, Tire2D, Tire3D + +PREFAB_BY_GEOMETRY = {Grid: ElasticGrid, Bar1D: ElasticBar, + Beam2D: ElasticBeam, Beam3D: ElasticBeam, + Tire: ElasticTire, Tire2D: ElasticTire, Tire3D: ElasticTire} diff --git a/common/sofa/prefabs/base.py b/common/sofa/prefabs/base.py new file mode 100644 index 0000000..1603d4a --- /dev/null +++ b/common/sofa/prefabs/base.py @@ -0,0 +1,62 @@ +"""Prefabs a study configures: they take the case blocks they need and build themselves.""" + +import Sofa +import Sofa.Core + + +def params(config, *exclude): + """Converts a set of parameters in a configuration dictionary to SOFA Data/parameterization. + Every key is passed on except 'type' and excluded ones marked by the user through *exclude.""" + skip = ('type',) + exclude + return {k: v for k, v in config.items() if k not in skip} + + +class ScenePrefab(Sofa.Prefab): + """Base SOFA prefab: holds `configs` and offers adders a subclass calls to build itself from it.""" + + def __init__(self, *args, configs=None, geometry=None, **kwargs): + # object.__setattr__: the C++ object Sofa's own __setattr__ needs doesn't exist yet, and + # Sofa.Prefab.__init__ below calls init(), which reads configs/geometry + object.__setattr__(self, 'configs', configs or {}) + object.__setattr__(self, 'geometry', geometry) + Sofa.Prefab.__init__(self, *args, **kwargs) + + def mechanical_node(self): + raise NotImplementedError + + def add_force_field(self, node, config, material, template): + """The *FEMForceField component under test, its Data merged from `config` and `material`""" + node.addObject(config['type'], name='fem', template=template, topology='@topology', + **params(material), **params(config)) + + def add_solvers(self, node, config): + """Newton (if present) + linear solver (+ preconditioner) + integration scheme.""" + + # Newton Solver + newton_config = config.get('newton') + if newton_config is not None: + node.addObject(newton_config['type'], name='newton', **params(newton_config)) + + # Linear Solver & Preconditioner (optional) + linear_config = config['linearSolver'] + precond_config = linear_config.get('preconditioner') + if precond_config is None: + node.addObject(linear_config['type'], name='linearSolver', **params(linear_config)) + else: + precond_system_config = precond_config['system'] + node.addObject(precond_system_config['type'], name='precondSystem', **params(precond_system_config)) + node.addObject(precond_config['type'], name='precond', linearSystem='@precondSystem', + **params(precond_config, 'system')) + + linear_system_config = linear_config['system'] + node.addObject(linear_system_config['type'], name='solverSystem', + preconditionerSystem='@precondSystem', **params(linear_system_config)) + node.addObject(linear_config['type'], name='linearSolver', linearSystem='@solverSystem', + preconditioner='@precond', **params(linear_config, 'system', 'preconditioner')) + + # ODE Integration Scheme + integration_scheme = params(config['integration']) + integration_scheme['linearSolver'] = '@linearSolver' + if newton_config is not None: + integration_scheme['newtonSolver'] = '@newton' + node.addObject(config['integration']['type'], name='ode', **integration_scheme) diff --git a/common/sofa/prefabs/grid.py b/common/sofa/prefabs/grid.py new file mode 100644 index 0000000..8af2906 --- /dev/null +++ b/common/sofa/prefabs/grid.py @@ -0,0 +1,89 @@ +"""SOFA prefabs: an elastic rectilinear structured grid in 1 to 3 dimensions, as a bar or a beam.""" + +from .base import ScenePrefab, params +from ..conventions import VEC_DIM, ELEMENTS +from ...geometry import Grid + +# Element kinds this prefab meshes, by the topological dimension of the geometry. +ELEMENT_KINDS = {1: ('edge',), 2: ('quad', 'tri'), 3: ('hexa', 'tet')} + + +def validate_parameters(name, config): + """Required ElasticGrid parameters, and an element kind its geometry can be meshed with.""" + required = ['geometry', 'grid_resolution', 'element'] + missing = [p for p in required if p not in config] + if missing: + raise ValueError(f"{name}: missing required parameters {missing}") + geometry = config['geometry'] + if not isinstance(geometry, Grid): + raise ValueError(f"{name}: geometry must be a Grid, " + f"got {type(geometry).__name__}") + kinds = ELEMENT_KINDS[geometry.dim] + if config['element'] not in kinds: + raise ValueError(f"{name}: element on a {geometry.dim}D {type(geometry).__name__} must be " + f"one of {list(kinds)}, got {config['element']!r}") + + +class ElasticGrid(ScenePrefab): + """Prefab for an elastic rectilinear structured grid model in SOFA.""" + + # The mechanical node's name in the scene graph; a named grid calls it what it is. + BODY = 'Body' + + prefabParameters = [ + {'name': 'grid_resolution', 'type': 'Vec3d', 'help': 'grid points per axis [nx, ny, nz]'}, + {'name': 'element', 'type': 'string', 'help': 'element kind (edge/tri/quad/tet/hexa)'}, + ] + + def __init__(self, *args, **kwargs): + validate_parameters(type(self).__name__, kwargs) + super().__init__(*args, **kwargs) + + def init(self): + dim = self.geometry.dim + extents = list(self.geometry.parameters) + [0.0] * (3 - dim) + VecType = VEC_DIM[self.geometry.spatial_dimensions] + element_kind = ELEMENTS[self.element.value] + grid_resolution = list(self.grid_resolution.value)[:dim] + [1] * (3 - dim) + + # Grid Topology Node + with self.addChild('Grid') as grid_node: + grid_node.addObject('RegularGridTopology', name='grid', + nx=int(grid_resolution[0]), ny=int(grid_resolution[1]), + nz=int(grid_resolution[2]), min=[0.0, 0.0, 0.0], max=extents) + + # Node containing the body's components + with self.addChild(self.BODY) as body: + self.body = body + + # Topology + if element_kind.grid_mapping is None: + body.addObject(element_kind.container, name='topology', position='@../Grid/grid.position', + **{element_kind.data_name: f'@../Grid/grid.{element_kind.data_name}'}) + else: + body.addObject(element_kind.container, name='topology', position='@../Grid/grid.position') + body.addObject(element_kind.grid_mapping, input='@../Grid/grid', output='@topology', + **params(self.configs.get('gridMapping', {}))) + body.addObject(element_kind.container.replace('Container', 'Modifier')) + # DOFs + body.addObject('MechanicalObject', name='dofs', template=VecType) + # The force field under test. + self.add_force_field(body, self.configs['forceField'], self.configs['material'], + template=f"{VecType},{element_kind.cpp}") + # ODE & Linear Solvers + self.add_solvers(body, self.configs['solvers']) + + def mechanical_node(self): + return self.body + + +class ElasticBar(ElasticGrid): + """Prefab for an elastic bar model in SOFA.""" + + BODY = 'Bar' + + +class ElasticBeam(ElasticGrid): + """Prefab for an elastic beam model in SOFA.""" + + BODY = 'Beam' diff --git a/common/sofa/prefabs/tire.py b/common/sofa/prefabs/tire.py new file mode 100644 index 0000000..0afb488 --- /dev/null +++ b/common/sofa/prefabs/tire.py @@ -0,0 +1,186 @@ +"""SOFA prefab: an elastic tire -- an annulus, extruded along z when it is 3D.""" + +import itertools + +import numpy as np + +from .base import ScenePrefab +from ..conventions import VEC_DIM, ELEMENTS +from ...geometry import Tire + +# Element kinds this prefab meshes, by the topological dimension of the geometry. +ELEMENT_KINDS = {2: ('quad', 'tri'), 3: ('hexa', 'tet')} + +# The two triangles one quadrilateral splits into, on one diagonal and on the other. +_TRIANGLE_PATTERNS = (((1, 2, 0), (3, 0, 2)), + ((0, 1, 3), (2, 3, 1))) + +# Swapping the two nodes of every hexahedron edge along one axis, as the permutation it applies to +# the eight corner slots. sofa::geometry::Hexahedron::{x,y,z}Edges. +_EDGE_SWAPS = ((1, 0, 3, 2, 5, 4, 7, 6), + (3, 2, 1, 0, 7, 6, 5, 4), + (4, 5, 6, 7, 0, 1, 2, 3)) + +# The six tetrahedra one hexahedron splits into, before and after an odd number of those swaps. +_TETRA_PATTERNS = { + False: ((0, 5, 1, 6), (0, 1, 3, 6), (1, 3, 6, 2), (6, 3, 0, 7), (6, 7, 0, 5), (7, 5, 4, 0)), + True: ((0, 5, 6, 1), (0, 1, 6, 3), (1, 3, 2, 6), (6, 3, 7, 0), (6, 7, 5, 0), (7, 5, 0, 4)), +} + + +def validate_parameters(config): + """Required ElasticTire parameters, and an element kind its geometry can be meshed with.""" + required = ['geometry', 'grid_resolution', 'element'] + missing = [p for p in required if p not in config] + if missing: + raise ValueError(f"ElasticTire: missing required parameters {missing}") + geometry = config['geometry'] + if not isinstance(geometry, Tire): + raise ValueError(f"ElasticTire: geometry must be a Tire2D or a Tire3D, " + f"got {type(geometry).__name__}") + kinds = ELEMENT_KINDS[geometry.dim] + if config['element'] not in kinds: + raise ValueError(f"ElasticTire: element on a {type(geometry).__name__} must be one of " + f"{list(kinds)}, got {config['element']!r}") + + +def annulus(geometry, cells): + """The annulus in the xy-plane: its nodes, and the quadrilaterals over them. + + `cells` counts cells radially and circumferentially. The ring closes on itself, so the + circumferential direction holds exactly that many nodes, not one more. + """ + radial, circumferential = cells + + radii = np.linspace(geometry.inner_radius, geometry.outer_radius, radial + 1) + angles = np.linspace(0.0, 2.0 * np.pi, circumferential, endpoint=False) + + r, theta = np.meshgrid(radii, angles, indexing='ij') + nodes = np.stack([r * np.cos(theta), r * np.sin(theta), np.zeros_like(r)], + axis=-1).reshape(-1, 3) + + node = np.arange(nodes.shape[0]).reshape(r.shape) + i, j = np.mgrid[0:radial, 0:circumferential] + ring = (j + 1) % circumferential # the last cell of the ring closes back onto j = 0 + + # (radial, circumferential) turns counter-clockwise, the winding GridTopology gives a cell. + corners = (node[i, j], node[i + 1, j], node[i + 1, ring], node[i, ring]) + return nodes, np.stack(corners, axis=-1).reshape(-1, 4) + + +def extrude(nodes, quads, length, axial): + """The annulus swept along z into `axial` layers: the swept nodes, and the hexahedra over them. + + Each layer repeats the annulus, so a hexahedron is its quadrilateral followed by the same one + a layer up -- again the winding GridTopology gives a cell. + """ + layer = nodes.shape[0] + swept = np.tile(nodes, (axial + 1, 1)) + swept[:, 2] = np.repeat(np.linspace(0.0, length, axial + 1), layer) + + base = np.arange(axial)[:, None, None] * layer + hexahedra = np.concatenate([quads[None] + base, quads[None] + base + layer], axis=-1) + return swept, hexahedra.reshape(-1, 8) + + +def triangulate(quads, cells): + """Two triangles per quadrilateral, split the way Quad2TriangleTopologicalMapping splits one. + + Neighbours share a whole edge, so any choice of diagonal conforms; alternating it cell by cell + only keeps the mesh from leaning the same way everywhere. + """ + i, j = np.mgrid[0:cells[0], 0:cells[1]] + diagonal = ((i ^ j) & 1).reshape(-1) + + triangles = np.empty((quads.shape[0], 2, 3), dtype=np.int64) + for choice, pattern in enumerate(_TRIANGLE_PATTERNS): + selected = diagonal == choice + triangles[selected] = quads[selected][:, pattern] + + return triangles.reshape(-1, 3) + + +def tetrahedralize(hexahedra, cells): + """Six tetrahedra per hexahedron, split the way Hexa2TetraTopologicalMapping splits one. + + The node labelling alternates with the parity of the cell's index, so neighbours agree on the + diagonal of every face they share. Across the seam that pairs the last circumferential cell + with the first, which only agree when there is an even number of them. + """ + radial, circumferential, axial = cells + if circumferential % 2: + raise ValueError(f"ElasticTire: an odd number of circumferential cells " + f"({circumferential}) leaves the tetrahedra of the closing cell " + f"disagreeing with the first on their shared diagonal") + + # Layer by layer, as `extrude` stacks them. + k, i, j = np.mgrid[0:axial, 0:radial, 0:circumferential] + parities = np.stack([i % 2 == 0, j % 2 == 1, k % 2 == 1], axis=-1).reshape(-1, 3) + + tetrahedra = np.empty((hexahedra.shape[0], 6, 4), dtype=np.int64) + for signature in itertools.product((False, True), repeat=3): + slots, swapped = list(range(8)), False + for axis, swap in enumerate(signature): + if swap: + slots = [slots[p] for p in _EDGE_SWAPS[axis]] + swapped = not swapped + pattern = [[slots[corner] for corner in tetra] for tetra in _TETRA_PATTERNS[swapped]] + + selected = np.all(parities == signature, axis=-1) + tetrahedra[selected] = hexahedra[selected][:, pattern] + + return tetrahedra.reshape(-1, 4) + + +def mesh(geometry, cells, element): + """The nodes of `geometry` and its elements of the given kind.""" + nodes, elements = annulus(geometry, cells[:2]) + if geometry.dim == 3: + nodes, elements = extrude(nodes, elements, geometry.extrusion_length, cells[2]) + + if element == 'tri': + elements = triangulate(elements, cells) + elif element == 'tet': + elements = tetrahedralize(elements, cells) + return nodes, elements + + +class ElasticTire(ScenePrefab): + """Prefab for an elastic tire model in SOFA.""" + + prefabParameters = [ + {'name': 'grid_resolution', 'type': 'Vec3d', + 'help': 'nodes per axis [radial, circumferential] and [, axial] when extruded'}, + {'name': 'element', 'type': 'string', 'help': 'element kind (quad/tri/hexa/tet)'}, + ] + + def __init__(self, *args, **kwargs): + validate_parameters(kwargs) + super().__init__(*args, **kwargs) + + def init(self): + VecType = VEC_DIM[self.geometry.spatial_dimensions] + element_kind = ELEMENTS[self.element.value] + # A study states nodes per axis, the mesh is laid out in cells. + cells = [int(count) - 1 for count in + list(self.grid_resolution.value)[:self.geometry.dim]] + + nodes, elements = mesh(self.geometry, cells, self.element.value) + + # Node containing Tire components + with self.addChild('Tire') as tire: + self.tire = tire + + # Topology: meshed directly, so no grid and no topological mapping to reach it. + tire.addObject(element_kind.container, name='topology', position=nodes.tolist(), + **{element_kind.data_name: elements.tolist()}) + # DOFs + tire.addObject('MechanicalObject', name='dofs', template=VecType) + # The force field under test. + self.add_force_field(tire, self.configs['forceField'], self.configs['material'], + template=f"{VecType},{element_kind.cpp}") + # ODE & Linear Solvers + self.add_solvers(tire, self.configs['solvers']) + + def mechanical_node(self): + return self.tire diff --git a/common/sofa/scene.py b/common/sofa/scene.py new file mode 100644 index 0000000..87de443 --- /dev/null +++ b/common/sofa/scene.py @@ -0,0 +1,29 @@ +"""Generic SOFA scene abstraction. Scenes in VnV build off of it.""" + +from abc import ABC, abstractmethod + + +class SofaScene(ABC): + """Abstract SOFA scene: loads its plugins, then builds the scene graph. + + A subclass appends the plugins it needs to `plugins`, and defines `body` to build the + scene graph. + """ + + plugins = [] # plugins every scene needs; a subclass appends its own on top + + def __init__(self): + # The list of plugins is kept as a copy per-instance. + self.plugins = list(self.plugins) + + def add_plugins(self, root): + """Add every plugin in `plugins` as one RequiredPlugin. Subclass appends to `self.plugins`.""" + root.addObject('RequiredPlugin', pluginName=self.plugins) + + @abstractmethod + def body(self, root): + """Build the scene graph.""" + + def build(self, root): + self.add_plugins(root) + return self.body(root) diff --git a/verification/__init__.py b/verification/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/verification/case.py b/verification/case.py new file mode 100644 index 0000000..eb06ae9 --- /dev/null +++ b/verification/case.py @@ -0,0 +1,209 @@ +"""The case as one object: parsed, validated, and resolved into a config dictionary that a run needs.""" + +import json +import pathlib + +from common.geometry import Grid, Bar1D, Beam2D, Beam3D, Tire, Tire2D, Tire3D +from common.sofa.conventions import ELEMENTS + +from . import manufactured +from .manufactured import ManufacturedProblem +from .manufactured.materials import LinearElastic +from .study.metrics import METRICS + +CASES_ROOT = pathlib.Path(__file__).parent / "cases" + +# case "geometry.type" -> geometry class, keyed by the class's own name. +GEOMETRIES = {cls.__name__: cls for cls in (Grid, Bar1D, Beam2D, Beam3D, + Tire, Tire2D, Tire3D)} + +# case "material.type" -> material class, keyed by the class's own name. +MATERIALS = {cls.__name__: cls for cls in (LinearElastic,)} + + +class _ManufacturedSolutionRegistry: + """(topological dim of the geometry, case "solution.type") -> ManufacturedSolution class. + + Looked up by naming convention: a case's solution type "Quadratic" at dim 2 resolves to + manufactured.Quadratic2D. The dimension is the PDE's own, not the space it is solved in. + """ + + def __getitem__(self, key): + dim, type_name = key + class_name = f"{type_name}{dim}D" + try: + return getattr(manufactured, class_name) + except AttributeError: + raise KeyError(f"no manufactured solution {class_name!r} for dim={dim}, " + f"solution type={type_name!r}") + + +MANUFACTURED_SOLUTIONS = _ManufacturedSolutionRegistry() + +# Required keys for all cases. A missing key is an omission, an unknown key is a typo. +REQUIRED = {"geometry", "element", "solution", "material", "forceField", + "errorQuadratureDegree", "sourceQuadratureDegree", "solvers", "mesh", + "asymptoticTolerance", "expectedOrder", "expectedOrderTolerance", "relativeNoiseFloor"} + +OPTIONAL = { + "tractionOn" # {region: directions} handed to the manufactured traction; the rest stay Dirichlet. + , "gridMapping" # Data for the mapping meshing the grid, where the element kind needs one. +} + +MESH_KEYS = {"cells", "levels", "refinementRatio"} + + + +def _validate(name, config): + """Every required key present, nothing unknown, and one expectation per metric. + + `solvers` and `forceField` pass their own keys through to SOFA Data verbatim, so this module + cannot know what's legal inside them and only checks their presence/shape, not their contents. + Every problem is collected before raising, so a case is corrected in one pass, not one key per run. + """ + problems = [] + + # Top-level keys: exact match against REQUIRED. + missing = REQUIRED - set(config) + if missing: + problems.append(f"missing {sorted(missing)}") + unknown = set(config) - REQUIRED - OPTIONAL + if unknown: + problems.append(f"unknown {sorted(unknown)}") + + # A grid-native element kind is meshed by no mapping, so its Data would go nowhere. + element = config.get("element") + if "gridMapping" in config and element in ELEMENTS and ELEMENTS[element].grid_mapping is None: + problems.append(f"gridMapping is stated but {element!r} is meshed straight from the grid") + + # mesh: exact match against MESH_KEYS. + if "mesh" in config and set(config["mesh"]) != MESH_KEYS: + problems.append(f"mesh must state exactly {sorted(MESH_KEYS)}, got {sorted(config['mesh'])}") + + # solution: just needs a type and an amplitude; the rest is that solution's own business. + solution_config = config.get("solution", {}) + absent = {"type", "amplitude"} - set(solution_config) + if absent: + problems.append(f"solution must state {sorted(absent)}") + + # material: must name a registered type. + material = config.get("material", {}) + if "type" not in material: + problems.append(f"material must state a type, one of {sorted(MATERIALS)}") + elif material["type"] not in MATERIALS: + problems.append(f"unknown material {material['type']!r}, expected one of {sorted(MATERIALS)}") + + # forceField: must be an object stating a type. + force_field = config.get("forceField") + if not isinstance(force_field, dict): + problems.append(f"forceField must be an object stating a type, got {force_field!r}") + force_field = {} + elif "type" not in force_field: + problems.append("forceField must state a type") + # forceField must not restate what the prefab already owns: name/template/topology follow + # from element and the mesh, so a case stating either could contradict what it already stated. + owned = {"name", "template", "topology"} & set(force_field) + if owned: + problems.append(f"forceField must not state {sorted(owned)}: the prefab sets them") + # forceField must not restate a constitutive parameter that belongs in material, so the + # manufactured source and the forceField under test can't drift apart. + clash = (set(force_field) & set(material)) - {"type"} + if clash: + problems.append(f"forceField and material both state {sorted(clash)}: state them in material") + + # expectedOrder: one entry per registered metric, exactly. + names = {metric.name for metric in METRICS} + if "expectedOrder" in config and set(config["expectedOrder"]) != names: + problems.append(f"expectedOrder must state exactly {sorted(names)}, " + f"got {sorted(config['expectedOrder'])}") + + if problems: + raise ValueError(f"{name}: " + "; ".join(problems)) + + +def _validate_boundary_conditions(name, config, geometry): + """The case's traction regions exist, and state one direction per dimension of the space.""" + traction_on = config.get("tractionOn", {}) + problems = [] + + unknown = set(traction_on) - set(geometry.region_names) + if unknown: + problems.append(f"tractionOn names {sorted(unknown)}, not regions of " + f"{type(geometry).__name__}: {sorted(geometry.region_names)}") + + # One direction per dimension of the mesh; the clamp pads the ones it is embedded in. + for region, directions in traction_on.items(): + if (len(directions) != geometry.dim + or any(direction not in (0, 1) for direction in directions)): + problems.append(f"tractionOn[{region!r}] must state {geometry.dim} directions, " + f"each 0 or 1, got {directions}") + + if problems: + raise ValueError(f"{name}: " + "; ".join(problems)) + + +class Case: + """One case: the study it specifies, with the geometry and the manufactured problem built.""" + + def __init__(self, path, config): + path = pathlib.Path(path) + + # Extract root path and keep stem as the name of the case. + stem = path.resolve().with_suffix("") + try: + self.name = stem.relative_to(CASES_ROOT.resolve()).as_posix() + except ValueError: + self.name = f"{path.parent.name}/{path.stem}" + _validate(self.name, config) + + # Geometry Class + geometry_config = dict(config["geometry"]) + self.geometry = GEOMETRIES[geometry_config.pop("type")](**geometry_config) + + _validate_boundary_conditions(self.name, config, self.geometry) + + # Keyed on the geometry's topological dimension. + self.traction_on = config.get("tractionOn", {}) + solution = MANUFACTURED_SOLUTIONS[(self.geometry.dim, config["solution"]["type"])]( + config["solution"], self.geometry, self.traction_on) + + # Nothing fixed anywhere leaves the rigid modes in the system. + if not solution.prescribe_displacement_on: + raise ValueError(f"{self.name}: tractionOn frees every direction of every region, " + f"which leaves the problem singular") + + # ManufacturedProblem Class + material = MATERIALS[config["material"]["type"]] + self.manufactured_problem = ManufacturedProblem( + solution, material(*material.get_material_parameters( + config["material"], self.geometry.spatial_dimensions)), + self.geometry.spatial_dimensions) + + self.element = config["element"] + self.material = config["material"] + self.force_field = config["forceField"] + self.solvers = config["solvers"] + self.mesh = config["mesh"] + self.error_quadrature_degree = config["errorQuadratureDegree"] + self.source_quadrature_degree = config["sourceQuadratureDegree"] + self.asymptotic_tolerance = config["asymptoticTolerance"] + self.expected_order = config["expectedOrder"] + self.expected_order_tolerance = config["expectedOrderTolerance"] + self.relative_noise_floor = config["relativeNoiseFloor"] + self.grid_mapping = config.get("gridMapping", {}) + + @classmethod + def load(cls, path): + with open(path) as f: + return cls(path, json.load(f)) + + @property + def equation(self): + """The manufactured solution in math form, for figures and reports.""" + return self.manufactured_problem.equation + + def levels(self): + """Cells per axis at each refinement level, coarsest first.""" + ratio = self.mesh["refinementRatio"] + return [[count * ratio ** level for count in self.mesh["cells"]] + for level in range(self.mesh["levels"])] diff --git a/verification/cases/LinearSmallStrainFEMForceField/Quadratic1D.json b/verification/cases/LinearSmallStrainFEMForceField/Quadratic1D.json new file mode 100644 index 0000000..2cfd0c2 --- /dev/null +++ b/verification/cases/LinearSmallStrainFEMForceField/Quadratic1D.json @@ -0,0 +1,24 @@ +{ + "geometry": { "type": "Bar1D", "length": 2.0 }, + "element": "edge", + "material": { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 }, + "forceField": { "type": "LinearSmallStrainFEMForceField" }, + "tractionOn": { "right": [1] }, + "solution": { "type": "Quadratic", "amplitude": 0.1 }, + "errorQuadratureDegree": 3, + "sourceQuadratureDegree": 2, + "solvers": { + "integration": { "type": "StaticSolver" }, + "newton": { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1, + "relativeInitialStoppingThreshold": 1e-6, + "relativeSuccessiveStoppingThreshold": 0, + "absoluteResidualStoppingThreshold": 0, + "printLog": false }, + "linearSolver": { "type": "SparseLDLSolver", "template": "CompressedRowSparseMatrixd" } + }, + "mesh": { "cells": [3], "levels": 7, "refinementRatio": 2 }, + "asymptoticTolerance": 0.10, + "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 }, + "expectedOrderTolerance": 0.15, + "relativeNoiseFloor": 1e-12 +} diff --git a/verification/cases/LinearSmallStrainFEMForceField/Quadratic2D_Quad.json b/verification/cases/LinearSmallStrainFEMForceField/Quadratic2D_Quad.json new file mode 100644 index 0000000..6b3824c --- /dev/null +++ b/verification/cases/LinearSmallStrainFEMForceField/Quadratic2D_Quad.json @@ -0,0 +1,24 @@ +{ + "geometry": { "type": "Beam2D", "length": 1.0, "width": 1.0 }, + "element": "quad", + "material": { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 }, + "forceField": { "type": "LinearSmallStrainFEMForceField" }, + "tractionOn": { "left": [0, 1], "bottom": [1, 0], "right": [1, 1], "top": [1, 1] }, + "solution": { "type": "Quadratic", "amplitude": 0.1 }, + "errorQuadratureDegree": 3, + "sourceQuadratureDegree": 3, + "solvers": { + "integration": { "type": "StaticSolver" }, + "newton": { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1, + "relativeInitialStoppingThreshold": 1e-6, + "relativeSuccessiveStoppingThreshold": 0, + "absoluteResidualStoppingThreshold": 0, + "printLog": false }, + "linearSolver": { "type": "SparseLDLSolver", "template": "CompressedRowSparseMatrixd" } + }, + "mesh": { "cells": [3, 3], "levels": 6, "refinementRatio": 2 }, + "asymptoticTolerance": 0.10, + "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 }, + "expectedOrderTolerance": 0.15, + "relativeNoiseFloor": 1e-12 +} diff --git a/verification/cases/LinearSmallStrainFEMForceField/Quadratic2D_Tri.json b/verification/cases/LinearSmallStrainFEMForceField/Quadratic2D_Tri.json new file mode 100644 index 0000000..a65c2c0 --- /dev/null +++ b/verification/cases/LinearSmallStrainFEMForceField/Quadratic2D_Tri.json @@ -0,0 +1,24 @@ +{ + "geometry": { "type": "Beam2D", "length": 1.0, "width": 1.0 }, + "element": "tri", + "material": { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 }, + "forceField": { "type": "LinearSmallStrainFEMForceField" }, + "tractionOn": { "left": [0, 1], "bottom": [1, 0], "right": [1, 1], "top": [1, 1] }, + "solution": { "type": "Quadratic", "amplitude": 0.1 }, + "errorQuadratureDegree": 2, + "sourceQuadratureDegree": 2, + "solvers": { + "integration": { "type": "StaticSolver" }, + "newton": { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1, + "relativeInitialStoppingThreshold": 1e-6, + "relativeSuccessiveStoppingThreshold": 0, + "absoluteResidualStoppingThreshold": 0, + "printLog": false }, + "linearSolver": { "type": "SparseLDLSolver", "template": "CompressedRowSparseMatrixd" } + }, + "mesh": { "cells": [3, 3], "levels": 6, "refinementRatio": 2 }, + "asymptoticTolerance": 0.10, + "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 }, + "expectedOrderTolerance": 0.15, + "relativeNoiseFloor": 1e-12 +} diff --git a/verification/cases/LinearSmallStrainFEMForceField/Quadratic3D_Hexa.json b/verification/cases/LinearSmallStrainFEMForceField/Quadratic3D_Hexa.json new file mode 100644 index 0000000..70df3b7 --- /dev/null +++ b/verification/cases/LinearSmallStrainFEMForceField/Quadratic3D_Hexa.json @@ -0,0 +1,31 @@ +{ + "geometry": { "type": "Beam3D", "length": 1.0, "width": 1.0, "height": 1.0 }, + "element": "hexa", + "material": { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 }, + "forceField": { "type": "LinearSmallStrainFEMForceField" }, + "tractionOn": { "left": [0, 1, 1], "bottom": [1, 0, 1], "front": [1, 1, 0], "right": [1, 1, 1], "top": [1, 1, 1], "back": [1, 1, 1] }, + "solution": { "type": "Quadratic", "amplitude": 0.1 }, + "errorQuadratureDegree": 3, + "sourceQuadratureDegree": 3, + "solvers": { + "integration": { "type": "StaticSolver" }, + "newton": { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1, + "relativeInitialStoppingThreshold": 1e-6, + "relativeSuccessiveStoppingThreshold": 0, + "absoluteResidualStoppingThreshold": 0, + "printLog": false }, + "linearSolver": { + "type": "PCGLinearSolver", "iterations": 1000, "tolerance": 1e-24, "printLog": false, + "system": { "type": "PreconditionedMatrixFreeSystem" }, + "preconditioner": { + "type": "SSORPreconditioner", + "system": { "type": "MatrixLinearSystem", "template": "CompressedRowSparseMatrixd" } + } + } + }, + "mesh": { "cells": [3, 3, 3], "levels": 5, "refinementRatio": 2 }, + "asymptoticTolerance": 0.10, + "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 }, + "expectedOrderTolerance": 0.15, + "relativeNoiseFloor": 1e-10 +} diff --git a/verification/cases/LinearSmallStrainFEMForceField/Quadratic3D_Tetra.json b/verification/cases/LinearSmallStrainFEMForceField/Quadratic3D_Tetra.json new file mode 100644 index 0000000..fa3ae31 --- /dev/null +++ b/verification/cases/LinearSmallStrainFEMForceField/Quadratic3D_Tetra.json @@ -0,0 +1,32 @@ +{ + "geometry": { "type": "Beam3D", "length": 1.0, "width": 1.0, "height": 1.0 }, + "element": "tet", + "material": { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 }, + "forceField": { "type": "LinearSmallStrainFEMForceField" }, + "tractionOn": { "left": [0, 1, 1], "bottom": [1, 0, 1], "front": [1, 1, 0], "right": [1, 1, 1], "top": [1, 1, 1], "back": [1, 1, 1] }, + "solution": { "type": "Quadratic", "amplitude": 0.1 }, + "errorQuadratureDegree": 2, + "sourceQuadratureDegree": 2, + "solvers": { + "integration": { "type": "StaticSolver" }, + "newton": { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1, + "relativeInitialStoppingThreshold": 1e-6, + "relativeSuccessiveStoppingThreshold": 0, + "absoluteResidualStoppingThreshold": 0, + "printLog": false }, + "linearSolver": { + "type": "PCGLinearSolver", "iterations": 1000, "tolerance": 1e-24, "printLog": false, + "system": { "type": "PreconditionedMatrixFreeSystem" }, + "preconditioner": { + "type": "SSORPreconditioner", + "system": { "type": "MatrixLinearSystem", "template": "CompressedRowSparseMatrixd" } + } + } + }, + "gridMapping": { "swapping": true }, + "mesh": { "cells": [4, 4, 4], "levels": 5, "refinementRatio": 2 }, + "asymptoticTolerance": 0.10, + "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 }, + "expectedOrderTolerance": 0.15, + "relativeNoiseFloor": 1e-10 +} diff --git a/verification/cases/LinearSmallStrainFEMForceField/Quartic2D_Quad.json b/verification/cases/LinearSmallStrainFEMForceField/Quartic2D_Quad.json new file mode 100644 index 0000000..f9a49bb --- /dev/null +++ b/verification/cases/LinearSmallStrainFEMForceField/Quartic2D_Quad.json @@ -0,0 +1,31 @@ +{ + "geometry": { "type": "Tire2D", "inner_radius": 0.25, "outer_radius": 0.5 }, + "element": "quad", + "material": { "type": "LinearElastic", "youngModulus": 1.0, "poissonRatio": 0.3 }, + "forceField": { "type": "LinearSmallStrainFEMForceField" }, + "tractionOn": { "outer": [1, 1] }, + "solution": { "type": "Quartic", "amplitude": 1.0 }, + "errorQuadratureDegree": 3, + "sourceQuadratureDegree": 3, + "solvers": { + "integration": { "type": "StaticSolver" }, + "newton": { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1, + "relativeInitialStoppingThreshold": 1e-6, + "relativeSuccessiveStoppingThreshold": 0, + "absoluteResidualStoppingThreshold": 0, + "printLog": false }, + "linearSolver": { + "type": "PCGLinearSolver", "iterations": 1000, "tolerance": 1e-24, "printLog": false, + "system": { "type": "PreconditionedMatrixFreeSystem" }, + "preconditioner": { + "type": "SSORPreconditioner", + "system": { "type": "MatrixLinearSystem", "template": "CompressedRowSparseMatrixd" } + } + } + }, + "mesh": { "cells": [2, 16], "levels": 5, "refinementRatio": 2 }, + "asymptoticTolerance": 0.10, + "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 }, + "expectedOrderTolerance": 0.15, + "relativeNoiseFloor": 1e-10 +} diff --git a/verification/cases/LinearSmallStrainFEMForceField/Quartic2D_Tri.json b/verification/cases/LinearSmallStrainFEMForceField/Quartic2D_Tri.json new file mode 100644 index 0000000..c161f73 --- /dev/null +++ b/verification/cases/LinearSmallStrainFEMForceField/Quartic2D_Tri.json @@ -0,0 +1,31 @@ +{ + "geometry": { "type": "Tire2D", "inner_radius": 0.25, "outer_radius": 0.5 }, + "element": "tri", + "material": { "type": "LinearElastic", "youngModulus": 1.0, "poissonRatio": 0.3 }, + "forceField": { "type": "LinearSmallStrainFEMForceField" }, + "tractionOn": { "outer": [1, 1] }, + "solution": { "type": "Quartic", "amplitude": 1.0 }, + "errorQuadratureDegree": 2, + "sourceQuadratureDegree": 2, + "solvers": { + "integration": { "type": "StaticSolver" }, + "newton": { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1, + "relativeInitialStoppingThreshold": 1e-6, + "relativeSuccessiveStoppingThreshold": 0, + "absoluteResidualStoppingThreshold": 0, + "printLog": false }, + "linearSolver": { + "type": "PCGLinearSolver", "iterations": 1000, "tolerance": 1e-24, "printLog": false, + "system": { "type": "PreconditionedMatrixFreeSystem" }, + "preconditioner": { + "type": "SSORPreconditioner", + "system": { "type": "MatrixLinearSystem", "template": "CompressedRowSparseMatrixd" } + } + } + }, + "mesh": { "cells": [2, 16], "levels": 5, "refinementRatio": 2 }, + "asymptoticTolerance": 0.10, + "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 }, + "expectedOrderTolerance": 0.15, + "relativeNoiseFloor": 1e-10 +} diff --git a/verification/cases/LinearSmallStrainFEMForceField/Quartic3D_Hexa.json b/verification/cases/LinearSmallStrainFEMForceField/Quartic3D_Hexa.json new file mode 100644 index 0000000..7bdde69 --- /dev/null +++ b/verification/cases/LinearSmallStrainFEMForceField/Quartic3D_Hexa.json @@ -0,0 +1,31 @@ +{ + "geometry": { "type": "Tire3D", "inner_radius": 0.25, "outer_radius": 0.5, "extrusion_length": 0.25 }, + "element": "hexa", + "material": { "type": "LinearElastic", "youngModulus": 1.0, "poissonRatio": 0.3 }, + "forceField": { "type": "LinearSmallStrainFEMForceField" }, + "tractionOn": { "front": [1, 1, 1], "back": [1, 1, 1] }, + "solution": { "type": "Quartic", "amplitude": 1.0 }, + "errorQuadratureDegree": 3, + "sourceQuadratureDegree": 3, + "solvers": { + "integration": { "type": "StaticSolver" }, + "newton": { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1, + "relativeInitialStoppingThreshold": 1e-6, + "relativeSuccessiveStoppingThreshold": 0, + "absoluteResidualStoppingThreshold": 0, + "printLog": false }, + "linearSolver": { + "type": "PCGLinearSolver", "iterations": 1000, "tolerance": 1e-24, "printLog": false, + "system": { "type": "PreconditionedMatrixFreeSystem" }, + "preconditioner": { + "type": "SSORPreconditioner", + "system": { "type": "MatrixLinearSystem", "template": "CompressedRowSparseMatrixd" } + } + } + }, + "mesh": { "cells": [2, 16, 2], "levels": 4, "refinementRatio": 2 }, + "asymptoticTolerance": 0.10, + "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 }, + "expectedOrderTolerance": 0.15, + "relativeNoiseFloor": 1e-10 +} diff --git a/verification/cases/LinearSmallStrainFEMForceField/Quartic3D_Tetra.json b/verification/cases/LinearSmallStrainFEMForceField/Quartic3D_Tetra.json new file mode 100644 index 0000000..d99fd10 --- /dev/null +++ b/verification/cases/LinearSmallStrainFEMForceField/Quartic3D_Tetra.json @@ -0,0 +1,32 @@ +{ + "geometry": { "type": "Tire3D", "inner_radius": 0.25, "outer_radius": 0.5, "extrusion_length": 0.25 }, + "element": "tet", + "material": { "type": "LinearElastic", "youngModulus": 1.0, "poissonRatio": 0.3 }, + "forceField": { "type": "LinearSmallStrainFEMForceField" }, + "tractionOn": { "front": [1, 1, 1], "back": [1, 1, 1] }, + "solution": { "type": "Quartic", "amplitude": 1.0 }, + "errorQuadratureDegree": 2, + "sourceQuadratureDegree": 2, + "solvers": { + "integration": { "type": "StaticSolver" }, + "newton": { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1, + "relativeInitialStoppingThreshold": 1e-6, + "relativeSuccessiveStoppingThreshold": 0, + "absoluteResidualStoppingThreshold": 0, + "printLog": false }, + "linearSolver": { + "type": "PCGLinearSolver", "iterations": 1000, "tolerance": 1e-24, "printLog": false, + "system": { "type": "PreconditionedMatrixFreeSystem" }, + "preconditioner": { + "type": "SSORPreconditioner", + "system": { "type": "MatrixLinearSystem", "template": "CompressedRowSparseMatrixd" } + } + } + }, + "gridMapping": { "swapping": true }, + "mesh": { "cells": [2, 16, 2], "levels": 4, "refinementRatio": 2 }, + "asymptoticTolerance": 0.10, + "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 }, + "expectedOrderTolerance": 0.15, + "relativeNoiseFloor": 1e-10 +} diff --git a/verification/cases/LinearSmallStrainFEMForceField/Trigonometric1D.json b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric1D.json new file mode 100644 index 0000000..827a2f2 --- /dev/null +++ b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric1D.json @@ -0,0 +1,24 @@ +{ + "geometry": { "type": "Bar1D", "length": 2.0 }, + "element": "edge", + "material": { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 }, + "forceField": { "type": "LinearSmallStrainFEMForceField" }, + "tractionOn": { "right": [1] }, + "solution": { "type": "Trigonometric", "amplitude": 1 }, + "errorQuadratureDegree": 3, + "sourceQuadratureDegree": 2, + "solvers": { + "integration": { "type": "StaticSolver" }, + "newton": { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1, + "relativeInitialStoppingThreshold": 1e-6, + "relativeSuccessiveStoppingThreshold": 0, + "absoluteResidualStoppingThreshold": 0, + "printLog": false }, + "linearSolver": { "type": "SparseLDLSolver", "template": "CompressedRowSparseMatrixd" } + }, + "mesh": { "cells": [3], "levels": 7, "refinementRatio": 2 }, + "asymptoticTolerance": 0.10, + "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 }, + "expectedOrderTolerance": 0.15, + "relativeNoiseFloor": 1e-12 +} diff --git a/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Quad.json b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Quad.json new file mode 100644 index 0000000..f608cf9 --- /dev/null +++ b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Quad.json @@ -0,0 +1,24 @@ +{ + "geometry": { "type": "Beam2D", "length": 1.0, "width": 1.0 }, + "element": "quad", + "material": { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 }, + "forceField": { "type": "LinearSmallStrainFEMForceField" }, + "tractionOn": { "left": [0, 1], "bottom": [1, 0], "right": [1, 1], "top": [1, 1] }, + "solution": { "type": "Trigonometric", "amplitude": 0.1 }, + "errorQuadratureDegree": 3, + "sourceQuadratureDegree": 3, + "solvers": { + "integration": { "type": "StaticSolver" }, + "newton": { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1, + "relativeInitialStoppingThreshold": 1e-6, + "relativeSuccessiveStoppingThreshold": 0, + "absoluteResidualStoppingThreshold": 0, + "printLog": false }, + "linearSolver": { "type": "SparseLDLSolver", "template": "CompressedRowSparseMatrixd" } + }, + "mesh": { "cells": [3, 3], "levels": 6, "refinementRatio": 2 }, + "asymptoticTolerance": 0.10, + "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 }, + "expectedOrderTolerance": 0.15, + "relativeNoiseFloor": 1e-12 +} diff --git a/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Quad_Dirichlet.json b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Quad_Dirichlet.json new file mode 100644 index 0000000..a57bd17 --- /dev/null +++ b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Quad_Dirichlet.json @@ -0,0 +1,23 @@ +{ + "geometry": { "type": "Beam2D", "length": 1.0, "width": 1.0 }, + "element": "quad", + "material": { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 }, + "forceField": { "type": "LinearSmallStrainFEMForceField" }, + "solution": { "type": "Trigonometric", "amplitude": 0.02 }, + "errorQuadratureDegree": 3, + "sourceQuadratureDegree": 3, + "solvers": { + "integration": { "type": "StaticSolver" }, + "newton": { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1, + "relativeInitialStoppingThreshold": 1e-6, + "relativeSuccessiveStoppingThreshold": 0, + "absoluteResidualStoppingThreshold": 0, + "printLog": false }, + "linearSolver": { "type": "SparseLDLSolver", "template": "CompressedRowSparseMatrixd" } + }, + "mesh": { "cells": [3, 3], "levels": 5, "refinementRatio": 2 }, + "asymptoticTolerance": 0.10, + "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 }, + "expectedOrderTolerance": 0.15, + "relativeNoiseFloor": 1e-12 +} diff --git a/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Tri.json b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Tri.json new file mode 100644 index 0000000..5024659 --- /dev/null +++ b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Tri.json @@ -0,0 +1,24 @@ +{ + "geometry": { "type": "Beam2D", "length": 1.0, "width": 1.0 }, + "element": "tri", + "material": { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 }, + "forceField": { "type": "LinearSmallStrainFEMForceField" }, + "tractionOn": { "left": [0, 1], "bottom": [1, 0], "right": [1, 1], "top": [1, 1] }, + "solution": { "type": "Trigonometric", "amplitude": 0.1 }, + "errorQuadratureDegree": 2, + "sourceQuadratureDegree": 2, + "solvers": { + "integration": { "type": "StaticSolver" }, + "newton": { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1, + "relativeInitialStoppingThreshold": 1e-6, + "relativeSuccessiveStoppingThreshold": 0, + "absoluteResidualStoppingThreshold": 0, + "printLog": false }, + "linearSolver": { "type": "SparseLDLSolver", "template": "CompressedRowSparseMatrixd" } + }, + "mesh": { "cells": [3, 3], "levels": 6, "refinementRatio": 2 }, + "asymptoticTolerance": 0.10, + "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 }, + "expectedOrderTolerance": 0.15, + "relativeNoiseFloor": 1e-12 +} diff --git a/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Tri_Dirichlet.json b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Tri_Dirichlet.json new file mode 100644 index 0000000..1efc304 --- /dev/null +++ b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Tri_Dirichlet.json @@ -0,0 +1,23 @@ +{ + "geometry": { "type": "Beam2D", "length": 1.0, "width": 1.0 }, + "element": "tri", + "material": { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 }, + "forceField": { "type": "LinearSmallStrainFEMForceField" }, + "solution": { "type": "Trigonometric", "amplitude": 0.02 }, + "errorQuadratureDegree": 2, + "sourceQuadratureDegree": 2, + "solvers": { + "integration": { "type": "StaticSolver" }, + "newton": { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1, + "relativeInitialStoppingThreshold": 1e-6, + "relativeSuccessiveStoppingThreshold": 0, + "absoluteResidualStoppingThreshold": 0, + "printLog": false }, + "linearSolver": { "type": "SparseLDLSolver", "template": "CompressedRowSparseMatrixd" } + }, + "mesh": { "cells": [3, 3], "levels": 5, "refinementRatio": 2 }, + "asymptoticTolerance": 0.10, + "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 }, + "expectedOrderTolerance": 0.15, + "relativeNoiseFloor": 1e-12 +} diff --git a/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Tri_PlaneStrain.json b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Tri_PlaneStrain.json new file mode 100644 index 0000000..c32ddef --- /dev/null +++ b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric2D_Tri_PlaneStrain.json @@ -0,0 +1,23 @@ +{ + "geometry": { "type": "Beam2D", "length": 1.0, "width": 1.0, "spatialDimensions": 3 }, + "element": "tri", + "material": { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 }, + "forceField": { "type": "LinearSmallStrainFEMForceField" }, + "solution": { "type": "Trigonometric", "amplitude": 0.1 }, + "errorQuadratureDegree": 2, + "sourceQuadratureDegree": 2, + "solvers": { + "integration": { "type": "StaticSolver" }, + "newton": { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1, + "relativeInitialStoppingThreshold": 1e-6, + "relativeSuccessiveStoppingThreshold": 0, + "absoluteResidualStoppingThreshold": 0, + "printLog": false }, + "linearSolver": { "type": "SparseLDLSolver", "template": "CompressedRowSparseMatrixd" } + }, + "mesh": { "cells": [3, 3], "levels": 6, "refinementRatio": 2 }, + "asymptoticTolerance": 0.10, + "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 }, + "expectedOrderTolerance": 0.15, + "relativeNoiseFloor": 1e-12 +} diff --git a/verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Hexa.json b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Hexa.json new file mode 100644 index 0000000..2e0c27f --- /dev/null +++ b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Hexa.json @@ -0,0 +1,31 @@ +{ + "geometry": { "type": "Beam3D", "length": 1.0, "width": 1.0, "height": 1.0 }, + "element": "hexa", + "material": { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 }, + "forceField": { "type": "LinearSmallStrainFEMForceField" }, + "tractionOn": { "left": [0, 1, 1], "bottom": [1, 0, 1], "front": [1, 1, 0], "right": [1, 1, 1], "top": [1, 1, 1], "back": [1, 1, 1] }, + "solution": { "type": "Trigonometric", "amplitude": 0.1 }, + "errorQuadratureDegree": 3, + "sourceQuadratureDegree": 3, + "solvers": { + "integration": { "type": "StaticSolver" }, + "newton": { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1, + "relativeInitialStoppingThreshold": 1e-6, + "relativeSuccessiveStoppingThreshold": 0, + "absoluteResidualStoppingThreshold": 0, + "printLog": false }, + "linearSolver": { + "type": "PCGLinearSolver", "iterations": 1000, "tolerance": 1e-24, "printLog": false, + "system": { "type": "PreconditionedMatrixFreeSystem" }, + "preconditioner": { + "type": "SSORPreconditioner", + "system": { "type": "MatrixLinearSystem", "template": "CompressedRowSparseMatrixd" } + } + } + }, + "mesh": { "cells": [3, 3, 3], "levels": 5, "refinementRatio": 2 }, + "asymptoticTolerance": 0.10, + "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 }, + "expectedOrderTolerance": 0.15, + "relativeNoiseFloor": 1e-10 +} diff --git a/verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Hexa_Dirichlet.json b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Hexa_Dirichlet.json new file mode 100644 index 0000000..f5c0ff8 --- /dev/null +++ b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Hexa_Dirichlet.json @@ -0,0 +1,30 @@ +{ + "geometry": { "type": "Beam3D", "length": 1.0, "width": 1.0, "height": 1.0 }, + "element": "hexa", + "material": { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 }, + "forceField": { "type": "LinearSmallStrainFEMForceField" }, + "solution": { "type": "Trigonometric", "amplitude": 0.02 }, + "errorQuadratureDegree": 3, + "sourceQuadratureDegree": 3, + "solvers": { + "integration": { "type": "StaticSolver" }, + "newton": { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1, + "relativeInitialStoppingThreshold": 1e-6, + "relativeSuccessiveStoppingThreshold": 0, + "absoluteResidualStoppingThreshold": 0, + "printLog": false }, + "linearSolver": { + "type": "PCGLinearSolver", "iterations": 1000, "tolerance": 1e-24, "printLog": false, + "system": { "type": "PreconditionedMatrixFreeSystem" }, + "preconditioner": { + "type": "SSORPreconditioner", + "system": { "type": "MatrixLinearSystem", "template": "CompressedRowSparseMatrixd" } + } + } + }, + "mesh": { "cells": [3, 3, 3], "levels": 5, "refinementRatio": 2 }, + "asymptoticTolerance": 0.10, + "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 }, + "expectedOrderTolerance": 0.15, + "relativeNoiseFloor": 1e-10 +} diff --git a/verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Tetra.json b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Tetra.json new file mode 100644 index 0000000..26c5075 --- /dev/null +++ b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Tetra.json @@ -0,0 +1,32 @@ +{ + "geometry": { "type": "Beam3D", "length": 1.0, "width": 1.0, "height": 1.0 }, + "element": "tet", + "material": { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 }, + "forceField": { "type": "LinearSmallStrainFEMForceField" }, + "tractionOn": { "left": [0, 1, 1], "bottom": [1, 0, 1], "front": [1, 1, 0], "right": [1, 1, 1], "top": [1, 1, 1], "back": [1, 1, 1] }, + "solution": { "type": "Trigonometric", "amplitude": 0.1 }, + "errorQuadratureDegree": 2, + "sourceQuadratureDegree": 2, + "solvers": { + "integration": { "type": "StaticSolver" }, + "newton": { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1, + "relativeInitialStoppingThreshold": 1e-6, + "relativeSuccessiveStoppingThreshold": 0, + "absoluteResidualStoppingThreshold": 0, + "printLog": false }, + "linearSolver": { + "type": "PCGLinearSolver", "iterations": 1000, "tolerance": 1e-24, "printLog": false, + "system": { "type": "PreconditionedMatrixFreeSystem" }, + "preconditioner": { + "type": "SSORPreconditioner", + "system": { "type": "MatrixLinearSystem", "template": "CompressedRowSparseMatrixd" } + } + } + }, + "gridMapping": { "swapping": true }, + "mesh": { "cells": [4, 4, 4], "levels": 5, "refinementRatio": 2 }, + "asymptoticTolerance": 0.10, + "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 }, + "expectedOrderTolerance": 0.15, + "relativeNoiseFloor": 1e-10 +} diff --git a/verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Tetra_Dirichlet.json b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Tetra_Dirichlet.json new file mode 100644 index 0000000..bb6b7ba --- /dev/null +++ b/verification/cases/LinearSmallStrainFEMForceField/Trigonometric3D_Tetra_Dirichlet.json @@ -0,0 +1,31 @@ +{ + "geometry": { "type": "Beam3D", "length": 1.0, "width": 1.0, "height": 1.0 }, + "element": "tet", + "material": { "type": "LinearElastic", "youngModulus": 1e6, "poissonRatio": 0.3 }, + "forceField": { "type": "LinearSmallStrainFEMForceField" }, + "solution": { "type": "Trigonometric", "amplitude": 0.02 }, + "errorQuadratureDegree": 2, + "sourceQuadratureDegree": 2, + "solvers": { + "integration": { "type": "StaticSolver" }, + "newton": { "type": "NewtonRaphsonSolver", "maxNbIterationsNewton": 1, + "relativeInitialStoppingThreshold": 1e-6, + "relativeSuccessiveStoppingThreshold": 0, + "absoluteResidualStoppingThreshold": 0, + "printLog": false }, + "linearSolver": { + "type": "PCGLinearSolver", "iterations": 1000, "tolerance": 1e-24, "printLog": false, + "system": { "type": "PreconditionedMatrixFreeSystem" }, + "preconditioner": { + "type": "SSORPreconditioner", + "system": { "type": "MatrixLinearSystem", "template": "CompressedRowSparseMatrixd" } + } + } + }, + "gridMapping": { "swapping": true }, + "mesh": { "cells": [4, 4, 4], "levels": 5, "refinementRatio": 2 }, + "asymptoticTolerance": 0.10, + "expectedOrder": { "L2": 2.0, "H1": 1.0, "EnergyDifference": 1.0 }, + "expectedOrderTolerance": 0.15, + "relativeNoiseFloor": 1e-10 +} diff --git a/verification/docs/IntroductionMMS.ipynb b/verification/docs/IntroductionMMS.ipynb new file mode 100644 index 0000000..8d5c668 --- /dev/null +++ b/verification/docs/IntroductionMMS.ipynb @@ -0,0 +1,889 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "01-md-title", + "metadata": {}, + "source": [ + "# Code Verification by the Method of Manufactured Solutions\n", + "\n", + "This notebook presents the Method of Manufactured Solutions (MMS) used for code verification and provides a walkthrough of the code that implements it.\n", + "\n", + "Verification and Validation answer two distinct questions [[2]](#ref2)[[6]](#ref6): \n", + "\n", + "- **Verification** Does the code solve its equations correctly, at the order of accuracy the discretisation claims?\n", + "- **Validation** Are these the right equations for the physical problem?\n", + "\n", + "In short, verification is *solving the equations right*, and validation is *solving the right equations* [[7]](#ref7).\n", + "\n", + "Verifying a code means comparing its solution against an exact one. The problem with classical exact solutions is that they exist\n", + "only for very simple geometries, loads and material laws, and even then they leave terms of the\n", + "equations untested [[2]](#ref2)[[4]](#ref4). A code can pass them and still be wrong. \n", + "MMS removes this limitation: it builds the problem around a solution we construct, so an exact\n", + "solution is always available." + ] + }, + { + "cell_type": "markdown", + "id": "02-md-mms-in-a-nutshell", + "metadata": {}, + "source": [ + "## The MMS in a nutshell\n", + "\n", + "Suppose we wish to solve the mathematical model $\\mathcal{L}(u) = 0$, where $\\mathcal{L}$ is the\n", + "differential operator and $u$ the unknown. The main idea is to choose a manufactured solution\n", + "$\\hat{u}$ first, then apply the governing equations to it: what remains, $s = \\mathcal{L}(\\hat{u})$,\n", + "is the source term. By construction, $\\hat{u}$ is the exact solution of the modified model\n", + "[[1]](#ref1)[[2]](#ref2)[[5]](#ref5)\n", + "\n", + "$$\\mathcal{L}(u) = s$$\n", + "\n", + "The boundary conditions are derived from $\\hat{u}$ as well. Three properties make this practical\n", + "[[5]](#ref5):\n", + "\n", + "- $s$ depends on the coordinates only, not on $u$. It enters the code as a known load, that is, as\n", + " input.\n", + "- Since $s$ is a known function, no discrete version of $\\mathcal{L}$ is applied to it. It is\n", + " only integrated into nodal loads, by quadrature.\n", + "- The left-hand side $\\mathcal{L}(u)$ is unchanged, so the discretisation under test is not\n", + " modified; only its input is.\n", + "\n", + "In short, the expectation is that solving the modified problem gives back $\\hat{u}$.\n", + "\n", + "In practice, the code does not solve $\\mathcal{L}(u) = s$ itself but a discrete version of it: a\n", + "finite element code solves its weak form on a mesh. Its solution $u_h$ is then only an\n", + "approximation of $\\hat{u}$, and the difference $e = u_h - \\hat{u}$ is the error. Because $\\hat{u}$\n", + "is known exactly, this error can be computed rather than estimated.\n", + "\n", + "### Measuring the error\n", + "\n", + "Numerical error has three sources [[3]](#ref3): the discretisation, finite-precision arithmetic\n", + "and the tolerance the solvers stop at. A coding mistake adds a fourth contribution, the one the\n", + "study is looking for:\n", + "\n", + "$$e = \\varepsilon_h + \\varepsilon_\\text{ro} + \\varepsilon_\\text{it} + \\varepsilon_\\text{bug}$$\n", + "\n", + "- $\\varepsilon_h$, the *discretisation error*, which comes from the mesh and vanishes as it is\n", + " refined;\n", + "- $\\varepsilon_\\text{ro}$, the *round-off error*, from finite-precision arithmetic;\n", + "- $\\varepsilon_\\text{it}$, the *iterative error*, from stopping the linear or Newton solver at a\n", + " finite tolerance;\n", + "- $\\varepsilon_\\text{bug}$, the error caused by *coding mistakes*.\n", + "\n", + "Under a correct implementation $\\varepsilon_\\text{bug} = 0$, and only $\\varepsilon_h$ should remain,\n", + "provided $\\varepsilon_\\text{ro}$ and $\\varepsilon_\\text{it}$ are kept far below it.\n", + "\n", + "The error is a field over $\\Omega$, so its size is measured with norms integrated over the domain:\n", + "\n", + "$$\\|e\\|_{L^2} = \\left(\\int_\\Omega |e|^2\\,d\\Omega\\right)^{1/2}, \\qquad\n", + " |e|_{H^1} = \\left(\\int_\\Omega \\|\\nabla e\\|^2\\,d\\Omega\\right)^{1/2}$$\n", + "\n", + "The $L^2$ norm measures the error itself, the $H^1$ semi-norm its gradient.\n", + "This suite also compares the strain energy SOFA computes for $u_h$ with the exact one, $E(\\hat{u})$, through $\\sqrt{2\\,|E(\\hat{u}) - E(u_h)|}$.\n", + "\n", + "Order verification assumes that the discretisation error can be expanded in powers of the mesh size\n", + "$h$ [[5]](#ref5):\n", + "\n", + "$$\\|\\varepsilon_h\\| = C\\,h^p + O(h^{p+1})$$\n", + "\n", + "$p$ is the formal order of accuracy of the discretisation. $C$ does not depend on $h$, but on the\n", + "higher derivatives of the manufactured solution [[5]](#ref5). As $h$ decreases, the higher-order\n", + "terms vanish faster than the first one, so on fine enough meshes the leading term dominates and\n", + "$\\|\\varepsilon_h\\| \\approx C\\,h^p$. Such meshes are in the *asymptotic regime*. On coarser ones, in the pre-asymptotic regime, the higher-order terms still weigh, and the observed order can differ from $p$ even for a correct code. For elements whose shape functions are polynomials of degree\n", + "$q$, $p = q + 1$ in the $L^2$ norm and $p = q$ in the $H^1$ semi-norm [[8]](#ref8).\n", + "\n", + "$\\varepsilon_\\text{ro}$ and $\\varepsilon_\\text{it}$ do not follow $h$: the iterative error is set\n", + "by the solver tolerance, and round-off grows as the mesh is refined and the system becomes larger\n", + "and worse conditioned. They set a noise floor, and an error below it must not be used to measure\n", + "the order: the ratio of two such errors reflects round-off and solver tolerance, not the\n", + "discretisation. Both must stay well below $\\varepsilon_h$ on the finest mesh [[5]](#ref5). A coding\n", + "mistake shows up as an error that does not vanish, or that vanishes at the wrong rate\n", + "[[4]](#ref4).\n", + "\n", + "The code is judged by the rate at which the error decreases, not by its size. A single error value\n", + "on a single mesh does not tell whether the code is correct: how small is small enough depends on the\n", + "solution, the units and the mesh, all of which enter $C$, while $p$ depends on the discretisation\n", + "alone. The mesh is therefore refined by a factor $r = h_\\text{coarse}/h_\\text{fine} > 1$, and the\n", + "ratio of the two errors gives the observed order of accuracy [[5]](#ref5):\n", + "\n", + "$$\\hat{p} = \\frac{\\ln\\!\\left(\\|e_{rh}\\| / \\|e_{h}\\|\\right)}{\\ln r}$$\n", + "\n", + "The test passes when $\\hat{p}$ matches $p$. Checking the order is a much stronger criterion than\n", + "checking that the error is small, or even that it decreases [[4]](#ref4).\n", + "\n", + "Everything under test is implemented in SOFA, in C++: the force fields\n", + "(`LinearSmallStrainFEMForceField`), the\n", + "solvers, and the components that turn a body force into nodal loads (`VectorSourceTerm`,\n", + "`FEMSourceTermIntegrator`). The suite does not modify any of them; they are the discrete\n", + "$\\mathcal{L}$ whose left-hand side stays unchanged.\n", + "\n", + "The Python scripts do everything around them:\n", + "\n", + "- derive the source term $s$ and the boundary conditions from $\\hat{u}$ and the material law;\n", + "- apply them to a SOFA scene as loads and constraints;\n", + "- solve the scene on a sequence of refined meshes, measure the error on each and read its order.\n", + "\n", + "### Further reading\n", + "\n", + "- Steinberg and Roache [[1]](#ref1): the original presentation of the method.\n", + "- Roache [[2]](#ref2): a short review of MMS for code verification.\n", + "- Roy [[3]](#ref3): a review of code and solution verification, with guidance on choosing manufactured solutions.\n", + "- Salari and Knupp [[4]](#ref4): the most complete treatment, including examples of coding mistakes detected by MSS\n", + "- Oberkampf and Roy [[5]](#ref5), Chapters 5 and 6: the textbook reference on code verification and MMS.\n", + "- Aycock *et al.* [[9]](#ref9): MMS applied to elastostatic solid mechanics in a finite element solver.\n" + ] + }, + { + "cell_type": "markdown", + "id": "03-md-example-problem", + "metadata": {}, + "source": [ + "## Example: a trigonometric manufactured solution on a 2D square\n", + "\n", + "The rest of this notebook walks through one verification study, building each piece of it in turn.\n", + "\n", + "### The governing equations\n", + "\n", + "Under small-strain deformations, the displacement $u$ of a linear elastic body in static equilibrium\n", + "and without external loads satisfies\n", + "\n", + "$$\\mathcal{L}(u) = -\\nabla\\cdot\\sigma(u) = 0, \\qquad\n", + " \\sigma = \\lambda\\,\\mathrm{tr}(\\varepsilon)\\,I + 2\\mu\\,\\varepsilon, \\qquad\n", + " \\varepsilon = \\tfrac{1}{2}\\left(\\nabla u + \\nabla u^T\\right)$$\n", + "\n", + "This is the equation solved by the `LinearSmallStrainFEMForceField` component in SOFA. $\\mu$ and\n", + "$\\lambda$ are the Lamé parameters, obtained from Young's modulus $E$ and Poisson's ratio $\\nu$.\n", + "\n", + "The body $\\Omega$ is a rectangle of length $L$ and width $W$, here a square with $L = W = 1$. Its\n", + "boundary is made of four sides: `left` ($x = 0$), `right` ($x = L$), `bottom` ($y = 0$) and `top`\n", + "($y = W$).\n", + "\n", + "### The manufactured solution\n", + "\n", + "We choose one sine period along each axis:\n", + "\n", + "$$\\hat{u}(x, y) = A \\begin{bmatrix} \\sin(kx)\\cos(ky) \\\\ \\cos(kx)\\sin(ky) \\end{bmatrix},\n", + " \\qquad k = \\frac{2\\pi}{L} = \\frac{2\\pi}{W}$$\n", + "\n", + "$\\hat{u}$ is not a physical solution and does not need to be one. It only has to be smooth, to leave\n", + "no derivative appearing in the equations identically zero, and to have derivatives small enough that\n", + "the asymptotic regime is reached on affordable meshes [[4]](#ref4)[[5]](#ref5). The amplitude $A$\n", + "sets the last of these.\n", + "\n", + "Its gradient\n", + "\n", + "$$\\hat{\\varepsilon} = \\nabla\\hat{u} = A k \\begin{bmatrix} c & -s \\\\ -s & c \\end{bmatrix},\n", + " \\qquad c = \\cos(kx)\\cos(ky), \\quad s = \\sin(kx)\\sin(ky)$$\n", + "\n", + "and the stress\n", + "\n", + "$$\\hat{\\sigma} = 2 A k \\begin{bmatrix} (\\lambda+\\mu)\\,c & -\\mu\\,s \\\\ -\\mu\\,s & (\\lambda+\\mu)\\,c \\end{bmatrix}$$\n", + "\n", + "Applying $\\mathcal{L}$ to $\\hat{u}$ gives the source term:\n", + "\n", + "$$s = \\mathcal{L}(\\hat{u}) = -\\nabla\\cdot\\hat{\\sigma} = 2k^2(\\lambda + 2\\mu)\\,\\hat{u}$$\n", + "\n", + "\n", + "### The boundary conditions\n", + "\n", + "The boundary conditions are taken from $\\hat{u}$ as well. Multiplying $\\mathcal{L}(u) = s$ by a test\n", + "function $v$, integrating over $\\Omega$ and integrating by parts gives the weak form\n", + "\n", + "$$\\int_\\Omega \\sigma : \\nabla v \\, d\\Omega - \\int_{\\partial\\Omega} t \\cdot v \\, d\\Gamma\n", + " = \\int_\\Omega s \\cdot v \\, d\\Omega, \\qquad t = \\sigma\\,n$$\n", + "\n", + "where $n$ is the outward normal and $t$ the traction. On each side, and in each direction, a\n", + "boundary condition prescribes one of two quantities:\n", + "\n", + "- the displacement, $u = \\hat{u}$ (Dirichlet). The test function vanishes there, and so does the\n", + " boundary term;\n", + "- the traction, $t = \\hat{t} = \\hat{\\sigma}\\,n$ (Neumann), which enters the boundary term as a known\n", + " load.\n", + "\n", + "Since $\\hat{u}$ is known everywhere, both conditions are consistent with it, and which one goes where\n", + "is a choice. Here, the displacement is prescribed on the whole boundary. It is not a zero Dirichlet condition, for instance on `left` it is $\\hat{u}(0, y) = (0,\\ A\\sin(ky))$.\n", + "\n", + "The modified problem is then\n", + "\n", + "$$-\\nabla\\cdot\\sigma(u) = 2k^2(\\lambda + 2\\mu)\\,\\hat{u} \\quad \\text{in } \\Omega, \\qquad\n", + " u = \\hat{u} \\quad \\text{on } \\partial\\Omega$$\n", + "\n", + "and $\\hat{u}$ is its exact solution. In this example $A = 0.02$, $E = 10^6$ and $\\nu = 0.3$." + ] + }, + { + "cell_type": "markdown", + "id": "04-md-problem-to-code", + "metadata": {}, + "source": [ + "### From the problem to the code\n", + "\n", + "Setting up an MMS study is systematic, so the problem above splits into separate concepts, each\n", + "with its own place in the repository:\n", + "\n", + "```\n", + "Verification-Validation/\n", + "├── common/ code shared by the suites\n", + "│ ├── geometry/ the domains Ω and their named boundary regions\n", + "│ └── sofa/\n", + "│ ├── prefabs/ how SOFA builds an elastic body on each geometry\n", + "│ ├── controllers.py SOFA controllers the scenes need, e.g. the prescribed displacement clamp\n", + "│ └── scene.py the base of every SOFA scene\n", + "└── verification/\n", + " ├── manufactured/ modules that implement manufactured solutions and material laws; combining them forms the manufactured problem to be solved\n", + " ├── scene/ the SOFA scene that loads a body with a manufactured problem\n", + " ├── study/ the error measurement and the convergence study\n", + " ├── case.py reads a case: the study's configuration\n", + " └── run.py runs a case\n", + "```\n", + "\n", + "Each of these is introduced below, where the walkthrough reaches it." + ] + }, + { + "cell_type": "markdown", + "id": "05-md-input-file", + "metadata": {}, + "source": [ + "### The input file\n", + "\n", + "A study is set up by an input file, a JSON file the suite reads into a dictionary. The walkthrough\n", + "builds it up as it goes, one key at a time. The problem above needs three keys:\n", + "\n", + "- `geometry` names the geometry and gives its size;\n", + "- `solution` names the manufactured solution and gives its amplitude $A$.\n", + "- `material` names the constitutive law and gives its parameters.\n", + "\n", + "Each `type` in the input file names a class." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "06-code-input-file", + "metadata": {}, + "outputs": [], + "source": [ + "input_file = {\n", + " \"geometry\": {\"type\": \"Beam2D\", \"length\": 1.0, \"width\": 1.0},\n", + " \"solution\": {\"type\": \"Trigonometric\", \"amplitude\": 0.02},\n", + " \"material\": {\"type\": \"LinearElastic\", \"youngModulus\": 1e6, \"poissonRatio\": 0.3},\n", + "}" + ] + }, + { + "cell_type": "markdown", + "id": "07-md-geometry", + "metadata": {}, + "source": [ + "### The manufactured problem\n", + "\n", + "The geometry is built first, because the solution depends on it." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "08-code-geometry", + "metadata": {}, + "outputs": [], + "source": [ + "import sys\n", + "\n", + "# Make the repository's packages importable from this folder.\n", + "sys.path.insert(0, \"../..\")\n", + "\n", + "from verification.case import GEOMETRIES\n", + "\n", + "geometry_parameters = dict(input_file[\"geometry\"])\n", + "geometry = GEOMETRIES[geometry_parameters.pop(\"type\")](**geometry_parameters)" + ] + }, + { + "cell_type": "markdown", + "id": "09-md-solution", + "metadata": {}, + "source": [ + "The solution's class is looked up by the geometry's dimension: `Trigonometric` in 2D is `Trigonometric2D` (`verification/manufactured/grid.py`). Solutions take the geometry as inputs because various parameters, such as $k$, might depend on the dimensions. Also, the displacement is prescribed on the geometry's regions." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10-code-solution", + "metadata": {}, + "outputs": [], + "source": [ + "from verification.case import MANUFACTURED_SOLUTIONS\n", + "\n", + "solution = MANUFACTURED_SOLUTIONS[(geometry.dim, input_file[\"solution\"][\"type\"])](input_file[\"solution\"], geometry)" + ] + }, + { + "cell_type": "markdown", + "id": "11-md-material", + "metadata": {}, + "source": [ + "The material converts the parameters of the input file into the ones its law is written with, which may depend on the number of spatial dimensions. `LinearElastic` (`verification/manufactured/materials.py`) converts $E$ and $\\nu$ to the Lamé parameters $\\mu$ and $\\lambda$; hyperelastic materials return their own." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12-code-material", + "metadata": {}, + "outputs": [], + "source": [ + "from verification.case import MATERIALS\n", + "\n", + "Material = MATERIALS[input_file[\"material\"][\"type\"]]\n", + "parameters = Material.get_material_parameters(input_file[\"material\"], geometry.spatial_dimensions)\n", + "material = Material(*parameters)" + ] + }, + { + "cell_type": "markdown", + "id": "13-md-problem", + "metadata": {}, + "source": [ + "The solution and the material are paired into the manufactured problem, which differentiates $\\hat{u}$ symbolically, with `sympy`, into the stress $\\hat{\\sigma}$ and the source term $s$. The problem is `ManufacturedProblem` (`verification/manufactured/problem.py`). It also compiles each expression into a numpy function evaluated at any array of points: `problem.u`, `problem.grad_u`, `problem.stress` and `problem.source`. These work as inputs to apply source terms such as body forces, prescribe displacement on boundaries and traction BCs. Additionally, the error is measured against them during the convergence analysis. These are the expressions it derived, with the numbers of the input file substituted." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14-code-problem", + "metadata": {}, + "outputs": [], + "source": [ + "import sympy as sp\n", + "from IPython.display import Math, display\n", + "\n", + "from verification.manufactured import ManufacturedProblem\n", + "\n", + "problem = ManufacturedProblem(solution, material, geometry.spatial_dimensions)\n", + "\n", + "\n", + "def rounded(expression):\n", + " \"\"\"The expression with its floating-point numbers cut to 6 significant digits, for display.\"\"\"\n", + " return expression.xreplace({f: sp.Float(f, 6) for f in expression.atoms(sp.Float)})\n", + "\n", + "\n", + "# The symbolic expressions the problem derived.\n", + "display(Math(r\"\\hat{u} = \" + sp.latex(rounded(problem.displacement_expression))),\n", + " Math(r\"\\hat{\\sigma} = \" + sp.latex(rounded(problem.stress_expression))),\n", + " Math(r\"s = \" + sp.latex(rounded(problem.source_expression))))" + ] + }, + { + "cell_type": "markdown", + "id": "15-md-plot-manufactured-solution", + "metadata": {}, + "source": [ + "The plot below uses `problem.u` to draw the exact $\\hat{u}$ over the square. On the right, the grid of the undeformed square is displaced by $\\hat{u}$: the boundary keeps its shape, but its points slide along it. The prescribed displacement is not zero on the boundary." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "16-code-plot-manufactured-solution", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "# û is small next to the square; magnify it so the deformation is visible.\n", + "MAGNIFICATION = 5\n", + "\n", + "# Evaluate the compiled û on a fine grid of points covering the square.\n", + "x, y = np.meshgrid(np.linspace(0, geometry.length, 201), np.linspace(0, geometry.width, 201))\n", + "u = problem.u(np.stack([x, y], axis=-1))\n", + "\n", + "fig, (magnitude, deformed) = plt.subplots(1, 2, figsize=(11, 4.8), layout=\"constrained\")\n", + "\n", + "# Left: the magnitude of û as colour, its direction as arrows on a coarser grid.\n", + "field = magnitude.pcolormesh(x, y, np.linalg.norm(u, axis=-1), shading=\"gouraud\")\n", + "fig.colorbar(field, ax=magnitude, label=r\"$|\\hat{u}|$\")\n", + "step = slice(None, None, 20)\n", + "magnitude.quiver(x[step, step], y[step, step], u[step, step, 0], u[step, step, 1], color=\"white\")\n", + "magnitude.set_title(r\"$\\hat{u}$: magnitude and direction\")\n", + "\n", + "# Right: each line of the undeformed grid, and the same line moved by the magnified û.\n", + "# The boundary lines are drawn thicker.\n", + "t = np.linspace(0, 1, 201)\n", + "for c in np.linspace(0, 1, 11):\n", + " for line in (np.stack([c * geometry.length + 0 * t, t * geometry.width], axis=-1),\n", + " np.stack([t * geometry.length, c * geometry.width + 0 * t], axis=-1)):\n", + " moved = line + MAGNIFICATION * problem.u(line)\n", + " deformed.plot(*line.T, color=\"0.8\", linewidth=0.6)\n", + " deformed.plot(*moved.T, color=\"black\", linewidth=1.6 if c in (0, 1) else 0.8)\n", + "deformed.set_title(rf\"Undeformed grid (grey) displaced by $\\hat{{u}}$, magnified ×{MAGNIFICATION}\")\n", + "\n", + "for ax in (magnitude, deformed):\n", + " ax.set_aspect(\"equal\")\n", + " ax.set_xlabel(\"$x$\")\n", + " ax.set_ylabel(\"$y$\")" + ] + }, + { + "cell_type": "markdown", + "id": "17-md-solving-keys", + "metadata": {}, + "source": [ + "### Solving on one mesh\n", + "\n", + "SOFA now solves the modified problem on a mesh. This needs four more keys for the `input_file`:\n", + "\n", + "- `forceField` is the SOFA component under test. It is the discrete $\\mathcal{L}$, the operator whose order of accuracy the study verifies;\n", + "- `element` is the kind of element the geometry is meshed with, here `quad`. SOFA builds the topology container from it;\n", + "- `solvers` sets how SOFA solves the discrete system;\n", + "- `sourceQuadratureDegree` is the degree of the quadrature that integrates the source term $s$ into nodal loads." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18-code-solving-keys", + "metadata": {}, + "outputs": [], + "source": [ + "input_file[\"forceField\"] = {\"type\": \"LinearSmallStrainFEMForceField\"}\n", + "input_file[\"element\"] = \"quad\"\n", + "input_file[\"solvers\"] = {\n", + " \"integration\": {\"type\": \"StaticSolver\"},\n", + " \"newton\": {\"type\": \"NewtonRaphsonSolver\", \"maxNbIterationsNewton\": 1,\n", + " \"relativeInitialStoppingThreshold\": 1e-6,\n", + " \"relativeSuccessiveStoppingThreshold\": 0,\n", + " \"absoluteResidualStoppingThreshold\": 0,\n", + " \"printLog\": False},\n", + " \"linearSolver\": {\"type\": \"SparseLDLSolver\", \"template\": \"CompressedRowSparseMatrixd\"},\n", + "}\n", + "input_file[\"sourceQuadratureDegree\"] = 3" + ] + }, + { + "cell_type": "markdown", + "id": "19-md-scene", + "metadata": {}, + "source": [ + "The scene is built in three steps:\n", + "\n", + "1. **The body.** A geometry knows nothing of SOFA; its prefab holds the instructions to build it in SOFA. `ElasticBeam` (`common/sofa/prefabs/grid.py`), the prefab of `Beam2D`, meshes the square with a regular grid of `quad` elements and adds the force field and the solvers.\n", + "2. **The source term.** `ApplyManufacturedSourceTerm` (`verification/scene/controllers.py`) evaluates `problem.source`, the $s$ derived above, at the nodes, and integrates it against the shape functions into nodal loads.\n", + "3. **The boundary conditions.** `RegionClamp` (`common/sofa/controllers.py`) adds one constraint per region where the solution prescribes the displacement, `problem.prescribe_displacement_on`, here all four sides. It moves their nodes to `problem.u`, the exact $\\hat{u}$, and fixes them there." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "20-code-scene", + "metadata": {}, + "outputs": [], + "source": [ + "import Sofa.Core\n", + "import Sofa.Simulation\n", + "\n", + "from common.sofa.conventions import VEC_DIM\n", + "from common.sofa.controllers import RegionClamp\n", + "from common.sofa.prefabs import PREFAB_BY_GEOMETRY\n", + "from verification.scene.controllers import ApplyManufacturedSourceTerm\n", + "from verification.scene.scene import MMSScene\n", + "\n", + "vec_type = VEC_DIM[geometry.spatial_dimensions]\n", + "\n", + "\n", + "def solve(cells):\n", + " \"\"\"Build the scene for a mesh of `cells` cells per side and solve it. Returns its root and body.\"\"\"\n", + " root = Sofa.Core.Node(\"root\")\n", + " root.addObject(\"RequiredPlugin\", pluginName=MMSScene.plugins)\n", + " root.addObject(\"DefaultAnimationLoop\")\n", + "\n", + " # 1.\n", + " Prefab = PREFAB_BY_GEOMETRY[type(geometry)]\n", + " prefab = root.addChild(Prefab(\n", + " name=Prefab.__name__, geometry=geometry, grid_resolution=[count + 1 for count in cells],\n", + " element=input_file[\"element\"],\n", + " configs={\"material\": input_file[\"material\"], \"forceField\": input_file[\"forceField\"],\n", + " \"solvers\": input_file[\"solvers\"]}))\n", + " body = prefab.mechanical_node()\n", + "\n", + " # 2.\n", + " body.addObject(ApplyManufacturedSourceTerm(\n", + " node=body, dofs=body.dofs, source=problem.source, vec_type=vec_type,\n", + " element=input_file[\"element\"], source_quadrature_degree=input_file[\"sourceQuadratureDegree\"],\n", + " name=\"bodyForceCtrl\"))\n", + "\n", + " # 3.\n", + " body.addObject(RegionClamp(\n", + " geometry=geometry, dofs=body.dofs, node=body,\n", + " prescribe_displacement_on=problem.prescribe_displacement_on,\n", + " vec_type=vec_type, displacement=problem.u, name=\"clampCtrl\"))\n", + "\n", + " # Initialise the scene, which runs the controllers, then solve it: one step of the static solver.\n", + " Sofa.Simulation.init(root)\n", + " Sofa.Simulation.animate(root, root.dt.value)\n", + " return root, body\n", + "\n", + "\n", + "# The mesh: 8 cells along each side of the square.\n", + "cells = [8, 8]\n", + "root, body = solve(cells)" + ] + }, + { + "cell_type": "markdown", + "id": "21-md-mms-scene", + "metadata": {}, + "source": [ + "Steps 1 to 3, building the scene for one mesh, are encapsulated in `MMSScene` (`verification/scene/scene.py`). It only builds the scene, initialization and solution are performed separately by the class that handles the mesh convergence study, introduced below. A side can also be given a traction instead of a displacement; the cases in `verification/cases/` show how." + ] + }, + { + "cell_type": "markdown", + "id": "22-md-plot-sofa-solution", + "metadata": {}, + "source": [ + "The plot draws SOFA's solution $u_h$, read off the nodes, over the exact $\\hat{u}$, both magnified as before." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "23-code-plot-sofa-solution", + "metadata": {}, + "outputs": [], + "source": [ + "from matplotlib.collections import PolyCollection\n", + "\n", + "# u_h at the nodes: where SOFA moved them, minus where they started.\n", + "nodes = body.dofs.rest_position.array()\n", + "u_h = body.dofs.position.array() - nodes\n", + "quads = body.topology.quads.array()\n", + "\n", + "fig, ax = plt.subplots(figsize=(5.5, 5.5), layout=\"constrained\")\n", + "\n", + "# The exact û, along the lines of the mesh.\n", + "for c in np.linspace(0, 1, cells[0] + 1):\n", + " for line in (np.stack([c * geometry.length + 0 * t, t * geometry.width], axis=-1),\n", + " np.stack([t * geometry.length, c * geometry.width + 0 * t], axis=-1)):\n", + " moved = line + MAGNIFICATION * problem.u(line)\n", + " ax.plot(*moved.T, color=\"0.6\", linewidth=0.8)\n", + "\n", + "# SOFA's solution u_h.\n", + "ax.add_collection(PolyCollection((nodes + MAGNIFICATION * u_h)[quads],\n", + " facecolor=\"none\", edgecolor=\"tab:red\", linewidth=1.2))\n", + "\n", + "ax.set_aspect(\"equal\")\n", + "ax.set_xlabel(\"$x$\")\n", + "ax.set_ylabel(\"$y$\")\n", + "ax.set_title(rf\"$u_h$ on the {cells[0]}×{cells[1]} mesh (red) and $\\hat{{u}}$ (grey), magnified ×{MAGNIFICATION}\")" + ] + }, + { + "cell_type": "markdown", + "id": "24-md-error", + "metadata": {}, + "source": [ + "### Measuring the error\n", + "\n", + "The norms $\\|e\\|_{L^2}$ and $|e|_{H^1}$ are integrated by quadrature over each element, with `MeshQuadrature` (`verification/study/quadrature.py`). At every quadrature point, $u_h$ and its gradient are interpolated from the nodes with the element's shape functions, whereas $\\hat{u}$ and $\\nabla\\hat{u}$ are evaluated exactly. One more key sets the quadrature: `errorQuadratureDegree`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "25-code-error", + "metadata": {}, + "outputs": [], + "source": [ + "from common.sofa.conventions import ELEMENTS\n", + "from verification.study.quadrature import MeshQuadrature\n", + "\n", + "input_file[\"errorQuadratureDegree\"] = 2\n", + "\n", + "\n", + "def error_norms(body):\n", + " \"\"\"‖e‖_L² and |e|_H¹ of the solution SOFA found on `body`'s mesh.\"\"\"\n", + " nodes = body.dofs.rest_position.array()\n", + " u_h = body.dofs.position.array() - nodes\n", + "\n", + " # The quadrature points and weights of every element of the mesh.\n", + " quadrature = MeshQuadrature(nodes, body.topology.quads.array(),\n", + " ELEMENTS[input_file[\"element\"]].cpp,\n", + " input_file[\"errorQuadratureDegree\"])\n", + "\n", + " # The error and its gradient at every quadrature point: u_h interpolated, û evaluated exactly.\n", + " e = quadrature.values(u_h) - quadrature.sample(problem.u)\n", + " grad_e = quadrature.grads(u_h) - quadrature.sample(problem.grad_u)\n", + "\n", + " # ‖e‖_L² = √(∫ |e|² dΩ) and |e|_H¹ = √(∫ ‖∇e‖² dΩ).\n", + " return (np.sqrt(quadrature.integrate(np.sum(e ** 2, axis=-1))),\n", + " np.sqrt(quadrature.integrate(np.sum(grad_e ** 2, axis=(-2, -1)))))\n", + "\n", + "\n", + "L2, H1 = error_norms(body)\n", + "print(f\"‖e‖_L2 = {L2:.3e}, |e|_H1 = {H1:.3e}\")" + ] + }, + { + "cell_type": "markdown", + "id": "26-md-metrics", + "metadata": {}, + "source": [ + "These two numbers alone do not say whether the code is correct. That is decided by how they decrease as the mesh is refined. The computation of the errors is encapsulated in `Measurement`, and each norm in a metric class: `L2`, `H1` and `EnergyDifference`. `METRICS` lists the metrics every mesh is measured with. All are in `verification/study/metrics.py`." + ] + }, + { + "cell_type": "markdown", + "id": "27-md-refinement", + "metadata": {}, + "source": [ + "### Refining the mesh\n", + "\n", + "The scene is now solved on a sequence of meshes, each $r$ times finer than the one before, and the\n", + "observed order $\\hat{p}$ is computed between each mesh and the one before it. The shape functions of `quad` elements\n", + "have degree $q = 1$, so they should give $p = 2$ in $L^2$ and $p = 1$ in $H^1$, and $p = 1$ for `EnergyDifference`. Three keys set this:\n", + "\n", + "- `mesh` gives the coarsest mesh's cells per side, the number of meshes, and $r$;\n", + "- `expectedOrder` gives $p$ for each metric;\n", + "- `expectedOrderTolerance` is how far $\\hat{p}$ on the finest pair of meshes may sit from $p$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "28-code-refinement", + "metadata": {}, + "outputs": [], + "source": [ + "# 3×3 cells, then 6×6, 12×12, 24×24 and 48×48.\n", + "input_file[\"mesh\"] = {\"cells\": [3, 3], \"levels\": 5, \"refinementRatio\": 2}\n", + "input_file[\"expectedOrder\"] = {\"L2\": 2.0, \"H1\": 1.0, \"EnergyDifference\": 1.0}\n", + "input_file[\"expectedOrderTolerance\"] = 0.15\n", + "\n", + "mesh = input_file[\"mesh\"]\n", + "spacings, norms = [], []\n", + "for level in range(mesh[\"levels\"]):\n", + " level_cells = [count * mesh[\"refinementRatio\"] ** level for count in mesh[\"cells\"]]\n", + " level_root, level_body = solve(level_cells)\n", + " spacings.append(geometry.length / level_cells[0])\n", + " norms.append(error_norms(level_body))\n", + " # Free this mesh's scene before building the next one.\n", + " Sofa.Simulation.unload(level_root)\n", + "\n", + "# One row per mesh, one column per norm: ‖e‖_L² and |e|_H¹.\n", + "norms = np.array(norms)\n", + "# p̂ = ln(‖e_rh‖ / ‖e_h‖) / ln r, between each mesh and the one before it.\n", + "orders = np.log(norms[:-1] / norms[1:]) / np.log(mesh[\"refinementRatio\"])\n", + "\n", + "print(f\"{'h':>8} {'L2':>10} {'p̂':>6} {'H1':>10} {'p̂':>6}\")\n", + "for level, h in enumerate(spacings):\n", + " rates = orders[level - 1] if level else (np.nan, np.nan)\n", + " print(f\"{h:8.4f} {norms[level, 0]:10.3e} {rates[0]:6.2f} {norms[level, 1]:10.3e} {rates[1]:6.2f}\")\n", + "\n", + "# The test: p̂ on the finest pair of meshes matches p, within the tolerance.\n", + "for name, order in zip((\"L2\", \"H1\"), orders[-1]):\n", + " expected = input_file[\"expectedOrder\"][name]\n", + " passed = abs(order - expected) <= input_file[\"expectedOrderTolerance\"]\n", + " print(f\"{name}: p̂ = {order:.2f}, p = {expected}: {'pass' if passed else 'fail'}\")" + ] + }, + { + "cell_type": "markdown", + "id": "29-md-convergence-study", + "metadata": {}, + "source": [ + "This loop is encapsulated in `ErrorConvergenceStudy` (`verification/study/study.py`). It judges the orders more strictly than the single check above:\n", + "\n", + "- a rate is discarded when either error is below the noise floor, `relativeNoiseFloor` times the size of the exact field. Below it, the round-off and iterative errors $\\varepsilon_\\text{ro}$ and $\\varepsilon_\\text{it}$ from the start dominate, and a rate measured there is meaningless;\n", + "- a rate is also discarded when the errors oscillate or diverge;\n", + "- the order is the one on the finest meshes, and is kept only once at least two consecutive rates agree within `asymptoticTolerance`: the meshes are then in the asymptotic regime, and rates from coarser meshes are discarded;\n", + "- the study fails if a Newton solve did not converge." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "30-code-convergence-study", + "metadata": {}, + "outputs": [], + "source": [ + "input_file[\"relativeNoiseFloor\"] = 1e-12\n", + "input_file[\"asymptoticTolerance\"] = 0.10" + ] + }, + { + "cell_type": "markdown", + "id": "31-md-full-input-file", + "metadata": {}, + "source": [ + "## Putting it together\n", + "\n", + "The input file is now complete. This is all a study needs, as the JSON file read by `run.py`:\n", + "\n", + "```json\n", + "{\n", + " \"geometry\": { \"type\": \"Beam2D\", \"length\": 1.0, \"width\": 1.0 },\n", + " \"element\": \"quad\",\n", + " \"material\": { \"type\": \"LinearElastic\", \"youngModulus\": 1e6, \"poissonRatio\": 0.3 },\n", + " \"forceField\": { \"type\": \"LinearSmallStrainFEMForceField\" },\n", + " \"solution\": { \"type\": \"Trigonometric\", \"amplitude\": 0.02 },\n", + " \"errorQuadratureDegree\": 3,\n", + " \"sourceQuadratureDegree\": 3,\n", + " \"solvers\": {\n", + " \"integration\": { \"type\": \"StaticSolver\" },\n", + " \"newton\": { \"type\": \"NewtonRaphsonSolver\", \"maxNbIterationsNewton\": 1,\n", + " \"relativeInitialStoppingThreshold\": 1e-6,\n", + " \"relativeSuccessiveStoppingThreshold\": 0,\n", + " \"absoluteResidualStoppingThreshold\": 0,\n", + " \"printLog\": false },\n", + " \"linearSolver\": { \"type\": \"SparseLDLSolver\", \"template\": \"CompressedRowSparseMatrixd\" }\n", + " },\n", + " \"mesh\": { \"cells\": [3, 3], \"levels\": 5, \"refinementRatio\": 2 },\n", + " \"asymptoticTolerance\": 0.10,\n", + " \"expectedOrder\": { \"L2\": 2.0, \"H1\": 1.0, \"EnergyDifference\": 1.0 },\n", + " \"expectedOrderTolerance\": 0.15,\n", + " \"relativeNoiseFloor\": 1e-12\n", + "}\n", + "```" + ] + }, + { + "cell_type": "markdown", + "id": "32-md-study", + "metadata": {}, + "source": [ + "Everything the walkthrough built by hand is done by two classes:\n", + "\n", + "- `Case` (`verification/case.py`) checks the input file, then builds the geometry, the solution, the material and the manufactured problem from it;\n", + "- `ErrorConvergenceStudy` builds the scene of every mesh with `MMSScene` and solves it, then computes the errors of that mesh's solution $u_h$ with `Measurement` and their norms with the metrics listed in `METRICS`, and finally judges the order of each metric. It prints one row per mesh with the rate of each metric, or why it was discarded, then the order each metric settled on next to the expected one.\n", + "\n", + "It also draws two figures: the error of each metric against $h$, its segments green where the rate was kept and red where it was discarded; and $|u_h|$ on each mesh, which approaches $|\\hat{u}|$ as the mesh is refined." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "33-code-study", + "metadata": {}, + "outputs": [], + "source": [ + "from verification.case import Case\n", + "from verification.study.study import ErrorConvergenceStudy\n", + "\n", + "# The path only names the case; the input file is passed as the dictionary built above.\n", + "case = Case(\"LinearSmallStrainFEMForceField/Trigonometric2D_Quad_Dirichlet.json\", input_file)\n", + "study = ErrorConvergenceStudy(case)\n", + "study.run()\n", + "print(\"\\nverified:\", study.verified)\n", + "\n", + "# The convergence plot the study drew while running: the error of each metric against h.\n", + "display(study.plots[\"convergence\"])\n", + "\n", + "# |u_h| on each mesh, coarsest to finest, drawn on the displaced mesh.\n", + "study.plot_fields()\n", + "display(study.plots[\"fields\"])" + ] + }, + { + "cell_type": "markdown", + "id": "34-md-running-a-case", + "metadata": {}, + "source": [ + "Saved as a JSON file under `verification/cases//`, the input file is a *case*, and the\n", + "runner does all of the above from the command line, writing the results to `verification/results/`:\n", + "\n", + "```bash\n", + "python verification/run.py LinearSmallStrainFEMForceField/Trigonometric2D_Quad_Dirichlet.json\n", + "python verification/run.py --all # every case, then one summary line per case\n", + "```\n", + "\n", + "To inspect one mesh of a case in the SOFA GUI, the runner also acts as a SOFA scene file:\n", + "\n", + "```bash\n", + "runSofa -l SofaPython3 verification/run.py --argv LinearSmallStrainFEMForceField/Trigonometric2D_Quad_Dirichlet.json --argv 0\n", + "```\n", + "\n", + "The first argument is the path to the .json case file and the second is the mesh refinement level to be built." + ] + }, + { + "cell_type": "markdown", + "id": "35-md-references", + "metadata": {}, + "source": [ + "\n", + "## References\n", + "\n", + "\n", + "**[1]** S. Steinberg and P. J. Roache, \"Symbolic manipulation and computational fluid dynamics,\"\n", + "*Journal of Computational Physics*, 57(2):251–284, 1985. See also P. J. Roache and\n", + "S. Steinberg, *AIAA Journal*, 22(10):1390–1394, 1984.\n", + "\n", + "\n", + "**[2]** P. J. Roache, \"Code verification by the method of manufactured solutions,\"\n", + "*Journal of Fluids Engineering*, 124(1):4–10, 2002.\n", + "[doi:10.1115/1.1436090](https://doi.org/10.1115/1.1436090)\n", + "\n", + "\n", + "**[3]** C. J. Roy, \"Review of code and solution verification procedures for computational\n", + "simulation,\" *Journal of Computational Physics*, 205(1):131–156, 2005.\n", + "[doi:10.1016/j.jcp.2004.10.036](https://doi.org/10.1016/j.jcp.2004.10.036)\n", + "\n", + "\n", + "**[4]** K. Salari and P. Knupp, *Code Verification by the Method of Manufactured Solutions*,\n", + "SAND2000-1444, Sandia National Laboratories, Albuquerque, NM, June 2000.\n", + "\n", + "\n", + "**[5]** W. L. Oberkampf and C. J. Roy, *Verification and Validation in Scientific Computing*,\n", + "Cambridge University Press, 2010. ISBN 978-0-521-11360-1. — *Chapter 5 (code verification and\n", + "order verification) and Chapter 6, §6.3 (MMS).*\n", + "\n", + "\n", + "**[6]** W. L. Oberkampf and T. G. Trucano, \"Verification and validation in computational fluid\n", + "dynamics,\" *Progress in Aerospace Sciences*, 38(3):209–272, 2002.\n", + "[doi:10.1016/S0376-0421(02)00005-2](https://doi.org/10.1016/S0376-0421%2802%2900005-2)\n", + "\n", + "\n", + "**[7]** B. W. Boehm, *Software Engineering Economics*, Prentice-Hall, Englewood Cliffs, NJ, 1981.\n", + "\n", + "\n", + "**[8]** P. G. Ciarlet, *The Finite Element Method for Elliptic Problems*, North-Holland, Amsterdam,\n", + "1978. Reprinted by SIAM, Classics in Applied Mathematics 40, 2002.\n", + "\n", + "\n", + "**[9]** K. I. Aycock, N. Rebelo and B. A. Craven, \"Method of manufactured solutions code\n", + "verification of elastostatic solid mechanics problems in a commercial finite element solver,\"\n", + "*Computers and Structures*, 229:106175, 2020.\n", + "[doi:10.1016/j.compstruc.2019.106175](https://doi.org/10.1016/j.compstruc.2019.106175)\n", + "\n", + "\n", + "**[10]** S. Étienne, A. Garon and D. Pelletier, \"Some manufactured solutions for verification of\n", + "fluid-structure interaction codes,\" *Computers and Structures*, 106–107:56–67, 2012.\n", + "[doi:10.1016/j.compstruc.2012.04.006](https://doi.org/10.1016/j.compstruc.2012.04.006)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/verification/docs/LinearSmallStrainFEMForceField.ipynb b/verification/docs/LinearSmallStrainFEMForceField.ipynb new file mode 100644 index 0000000..38fe22c --- /dev/null +++ b/verification/docs/LinearSmallStrainFEMForceField.ipynb @@ -0,0 +1,497 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "01-md-title", + "metadata": {}, + "source": [ + "# Verification of `LinearSmallStrainFEMForceField`\n", + "\n", + "This notebook verifies SOFA's `LinearSmallStrainFEMForceField` with the Method of Manufactured Solutions, presented in [IntroductionMMS.ipynb](IntroductionMMS.ipynb). It presents the governing PDEs and some of the cases used for verification, in `verification/cases/LinearSmallStrainFEMForceField/`." + ] + }, + { + "cell_type": "markdown", + "id": "02-md-material", + "metadata": {}, + "source": [ + "## Governing PDE for Linear Elastic Materials\n", + "\n", + "Under small-strain deformations, the displacement $\\mathbf{u}$ of a linear elastic body satisfies\n", + "\n", + "$$\\mathcal{L}(\\mathbf{u}) = -\\nabla \\cdot \\boldsymbol{\\sigma}(\\mathbf{u}) = 0$$\n", + "\n", + "This is the equation solved by the `LinearSmallStrainFEMForceField` component in SOFA. The body occupies a domain $\\Omega$ in $d = 1$, $2$ or $3$ dimensions, and $\\Gamma$ is its boundary.\n", + "\n", + "The material is defined by its strain energy density $\\psi$, the elastic energy stored per unit volume, and the stress is its derivative with respect to the displacement gradient. For a linear elastic material, $\\psi$ depends on the strain $\\boldsymbol{\\varepsilon}$:\n", + "\n", + "$$\\boldsymbol{\\varepsilon}(\\mathbf{u}) = \\frac{1}{2}\\left(\\nabla \\mathbf{u} + \\nabla \\mathbf{u}^T\\right), \\qquad\n", + " \\psi(\\boldsymbol{\\varepsilon}) = \\frac{\\lambda}{2}\\,\\mathrm{tr}(\\boldsymbol{\\varepsilon})^2 + \\mu\\,\\boldsymbol{\\varepsilon} : \\boldsymbol{\\varepsilon}$$\n", + "\n", + "$$\\boldsymbol{\\sigma}(\\mathbf{u}) = \\frac{\\partial \\psi}{\\partial \\nabla \\mathbf{u}} = \\lambda\\,\\mathrm{tr}(\\boldsymbol{\\varepsilon})\\,\\mathbf{I} + 2\\mu\\,\\boldsymbol{\\varepsilon}$$\n", + "\n", + "The Lamé parameters $\\lambda$ and $\\mu$ are computed from Young's modulus $E$ and the Poisson ratio $\\nu$:\n", + "\n", + "$$\\mu = \\frac{E}{2(1 + \\nu)}, \\qquad \\lambda = \\frac{E\\nu}{(1 + \\nu)\\,(1 - (d - 1)\\,\\nu)}$$\n", + "\n", + "### Manufactured Solution and Boundary Conditions\n", + "\n", + "For a manufactured solution $\\hat{\\mathbf{u}}$, with stress $\\hat{\\boldsymbol{\\sigma}} = \\boldsymbol{\\sigma}(\\hat{\\mathbf{u}})$, the source term is\n", + "\n", + "$$\\mathbf{s} = \\mathcal{L}(\\hat{\\mathbf{u}}) = -\\nabla \\cdot \\hat{\\boldsymbol{\\sigma}}$$\n", + "\n", + "The boundary conditions are taken from $\\hat{\\mathbf{u}}$ as well. The weak form shows which quantities a boundary condition can prescribe [[1]](#ref1). Multiplying $\\mathcal{L}(\\mathbf{u}) = \\mathbf{s}$ by a test function $\\mathbf{v}$, integrating over $\\Omega$ and applying the divergence theorem gives\n", + "\n", + "$$\\int_\\Omega \\boldsymbol{\\sigma} : \\boldsymbol{\\varepsilon}(\\mathbf{v})\\,d\\Omega - \\int_\\Gamma \\mathbf{t} \\cdot \\mathbf{v}\\,d\\Gamma\n", + " = \\int_\\Omega \\mathbf{s} \\cdot \\mathbf{v}\\,d\\Omega$$\n", + "\n", + "where $\\mathbf{t} = \\boldsymbol{\\sigma} \\cdot \\mathbf{n}$ is the traction and $\\mathbf{n}$ the outward normal of $\\Gamma$.\n", + "\n", + "At each point of $\\Gamma$ and in each direction, a boundary condition prescribes one of two quantities:\n", + "\n", + "- the displacement, $\\mathbf{u} = \\hat{\\mathbf{u}}$ (Dirichlet). The test function vanishes there, and so does the boundary term;\n", + "- the traction, $\\mathbf{t} = \\hat{\\mathbf{t}} = \\hat{\\boldsymbol{\\sigma}} \\cdot \\mathbf{n}$ (Neumann), which enters the boundary term as a known load.\n", + "\n", + "Since $\\hat{\\mathbf{u}}$ is known everywhere, both $\\hat{\\mathbf{u}}$ and $\\hat{\\mathbf{t}}$ are known on the whole of $\\Gamma$. Which one goes where is a choice, and a region of $\\Gamma$ can have its displacement prescribed in one direction and its traction in another; enough of the displacement must be prescribed to keep the body from translating and rotating freely.\n", + "\n", + "The modified problem is then\n", + "\n", + "$$-\\nabla \\cdot \\boldsymbol{\\sigma}(\\mathbf{u}) = \\mathbf{s} \\quad \\text{in } \\Omega, \\qquad\n", + " \\mathbf{u} = \\hat{\\mathbf{u}} \\quad \\text{on } \\Gamma_D, \\qquad\n", + " \\boldsymbol{\\sigma}(\\mathbf{u}) \\cdot \\mathbf{n} = \\hat{\\mathbf{t}} \\quad \\text{on } \\Gamma_N$$\n", + "\n", + "with $\\Gamma_D$ and $\\Gamma_N$ the parts of the boundary where the displacement and the traction are prescribed, and $\\hat{\\mathbf{u}}$ is its exact solution. Each case below picks $\\hat{\\mathbf{u}}$, $\\Omega$, $\\Gamma_D$ and $\\Gamma_N$." + ] + }, + { + "cell_type": "markdown", + "id": "03-md-setup", + "metadata": {}, + "source": [ + "## Running the Cases\n", + "\n", + "The cell below loads the suite's runner and prepares the figures. Run it once before any case: every case cell below then runs on its own, in any order." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "04-code-setup", + "metadata": {}, + "outputs": [], + "source": [ + "# Make the plugin's packages importable from verification/docs\n", + "import sys\n", + "sys.path.insert(0, \"../..\")\n", + "\n", + "# Hide the runner's printed report; figures still show\n", + "import os\n", + "sys.stdout = open(os.devnull, \"w\")\n", + "\n", + "# Load the suite's runner\n", + "from verification.run import run\n", + "\n", + "# For showing figures inline\n", + "%matplotlib inline\n", + "\n", + "# For drawing several convergence figures side by side\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "markdown", + "id": "05-md-1d", + "metadata": {}, + "source": [ + "## Edge Elements on a 1D Bar\n", + "\n", + "The bar starts at the origin and has length $L$ along $x$. The displacement $u$ is a scalar." + ] + }, + { + "cell_type": "markdown", + "id": "06-md-quadratic-1d", + "metadata": {}, + "source": [ + "### Quadratic Function\n", + "\n", + "$$\\hat{u} = A\\,x^2, \\qquad \\hat{\\sigma} = 2EAx, \\qquad s = -2EA$$\n", + "\n", + "- `left` ($x = 0$): $\\hat{u} = 0$\n", + "- `right` ($x = L$): $\\hat{t} = 2EAL$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "07-code-quadratic-1d", + "metadata": {}, + "outputs": [], + "source": [ + "# Case: 1D bar, quadratic manufactured solution\n", + "study = run(\"../cases/LinearSmallStrainFEMForceField/Quadratic1D.json\")\n", + "# u_h along the bar at each level, coarse to fine\n", + "study.plot_fields()\n", + "study.plots[\"convergence\"]" + ] + }, + { + "cell_type": "markdown", + "id": "08-md-trigonometric-1d", + "metadata": {}, + "source": [ + "### Trigonometric Function\n", + "\n", + "$$\\hat{u} = A\\sin(kx), \\qquad \\hat{\\sigma} = EAk\\cos(kx), \\qquad s = EAk^2\\sin(kx)$$\n", + "\n", + "with $k = \\frac{2\\pi}{L}$.\n", + "\n", + "- `left` ($x = 0$): $\\hat{u} = 0$\n", + "- `right` ($x = L$): $\\hat{t} = EAk$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "09-code-trigonometric-1d", + "metadata": {}, + "outputs": [], + "source": [ + "# Case: 1D bar, trigonometric manufactured solution\n", + "study = run(\"../cases/LinearSmallStrainFEMForceField/Trigonometric1D.json\")\n", + "# u_h along the bar at each level, coarse to fine\n", + "study.plot_fields()\n", + "study.plots[\"convergence\"]" + ] + }, + { + "cell_type": "markdown", + "id": "10-md-2d", + "metadata": {}, + "source": [ + "## Triangle and Quadrilateral Elements on a 2D Beam\n", + "\n", + "The beam starts at the origin, with length $L$ along $x$ and width $W$ along $y$." + ] + }, + { + "cell_type": "markdown", + "id": "11-md-quadratic-2d", + "metadata": {}, + "source": [ + "### Quadratic Function\n", + "\n", + "$$\n", + "\\hat{\\mathbf{u}} = A\\begin{pmatrix} x^2 + xy \\\\ y^2 + xy \\end{pmatrix}, \n", + "\\qquad \\mathbf{s} = -A \\begin{pmatrix} 3\\lambda + 5\\mu \\\\ 3\\lambda + 5\\mu \\end{pmatrix}, \n", + "\\qquad \\hat{\\boldsymbol{\\sigma}} = A\\begin{pmatrix} 3\\lambda(x + y) + 2\\mu(2x + y) & \\mu(x + y) \\\\ \\mu(x + y) & 3\\lambda(x + y) + 2\\mu(x + 2y) \\end{pmatrix}\n", + "$$\n", + "\n", + "- `left` ($x = 0$): $\\hat{u}_x = 0$, $\\hat{t}_y = -\\mu A y$\n", + "- `bottom` ($y = 0$): $\\hat{u}_y = 0$, $\\hat{t}_x = -\\mu A x$\n", + "- `right` ($x = L$): $\\hat{\\mathbf{t}} = A\\,\\big(3\\lambda(L + y) + 2\\mu(2L + y),\\ \\mu(L + y)\\big)$\n", + "- `top` ($y = W$): $\\hat{\\mathbf{t}} = A\\,\\big(\\mu(x + W),\\ 3\\lambda(x + W) + 2\\mu(x + 2W)\\big)$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "12-code-quadratic-2d", + "metadata": {}, + "outputs": [], + "source": [ + "# Case: 2D rectangle, quadratic manufactured solution, triangles (left) and quadrilaterals (right)\n", + "tri = run(\"../cases/LinearSmallStrainFEMForceField/Quadratic2D_Tri.json\")\n", + "quad = run(\"../cases/LinearSmallStrainFEMForceField/Quadratic2D_Quad.json\")\n", + "figure, (triangles, quadrilaterals) = plt.subplots(1, 2, figsize=(14.4, 4.8))\n", + "tri.plot(triangles)\n", + "quad.plot(quadrilaterals)\n", + "\n", + "# Deformed shape at each level, coarse to fine\n", + "tri.plot_fields()\n", + "quad.plot_fields()" + ] + }, + { + "cell_type": "markdown", + "id": "13-md-trigonometric-2d", + "metadata": {}, + "source": [ + "### Trigonometric Function\n", + "\n", + "$$\\hat{\\mathbf{u}} = A\\begin{pmatrix} \\sin(kx)\\cos(ky) \\\\ \\cos(kx)\\sin(ky) \\end{pmatrix}, \\qquad \\mathbf{s} = 2k^2(\\lambda + 2\\mu)\\,\\hat{\\mathbf{u}}, \\qquad \\hat{\\boldsymbol{\\sigma}} = 2Ak\\begin{pmatrix} (\\lambda + \\mu)\\cos(kx)\\cos(ky) & -\\mu\\sin(kx)\\sin(ky) \\\\ -\\mu\\sin(kx)\\sin(ky) & (\\lambda + \\mu)\\cos(kx)\\cos(ky) \\end{pmatrix}$$\n", + "\n", + "with $k = \\frac{2\\pi}{L}$, the same along $x$ and $y$ since the cases use a square, $L = W$.\n", + "\n", + "- `left` ($x = 0$): $\\hat{u}_x = 0$, $\\hat{t}_y = 0$\n", + "- `bottom` ($y = 0$): $\\hat{u}_y = 0$, $\\hat{t}_x = 0$\n", + "- `right` ($x = L$): $\\hat{\\mathbf{t}} = 2Ak(\\lambda + \\mu)\\,\\big(\\cos(ky),\\ 0\\big)$\n", + "- `top` ($y = L$): $\\hat{\\mathbf{t}} = 2Ak(\\lambda + \\mu)\\,\\big(0,\\ \\cos(kx)\\big)$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "14-code-trigonometric-2d", + "metadata": {}, + "outputs": [], + "source": [ + "# Case: 2D rectangle, trigonometric manufactured solution, triangles (left) and quadrilaterals (right)\n", + "tri = run(\"../cases/LinearSmallStrainFEMForceField/Trigonometric2D_Tri.json\")\n", + "quad = run(\"../cases/LinearSmallStrainFEMForceField/Trigonometric2D_Quad.json\")\n", + "figure, (triangles, quadrilaterals) = plt.subplots(1, 2, figsize=(14.4, 4.8))\n", + "tri.plot(triangles)\n", + "quad.plot(quadrilaterals)\n", + "\n", + "# Deformed shape at each level, coarse to fine\n", + "tri.plot_fields()\n", + "quad.plot_fields()" + ] + }, + { + "cell_type": "markdown", + "id": "15-md-trigonometric-2d-plane-strain", + "metadata": {}, + "source": [ + "#### Plane Strain\n", + "\n", + "The beam is solved in 3D space with $u_z = 0$. $\\lambda$ and $\\mu$ take their 3D values. The displacement is prescribed on all four sides:\n", + "\n", + "- `left` ($x = 0$) and `right` ($x = L$): $\\hat{\\mathbf{u}} = A\\,\\big(0,\\ \\sin(ky),\\ 0\\big)$\n", + "- `bottom` ($y = 0$) and `top` ($y = L$): $\\hat{\\mathbf{u}} = A\\,\\big(\\sin(kx),\\ 0,\\ 0\\big)$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "16-code-trigonometric-2d-plane-strain-tri", + "metadata": {}, + "outputs": [], + "source": [ + "# Case: 2D rectangle in plane strain, trigonometric manufactured solution, triangles\n", + "study = run(\"../cases/LinearSmallStrainFEMForceField/Trigonometric2D_Tri_PlaneStrain.json\")\n", + "# Deformed shape at each level, coarse to fine\n", + "study.plot_fields()\n", + "study.plots[\"convergence\"]" + ] + }, + { + "cell_type": "markdown", + "id": "17-md-2d-tire", + "metadata": {}, + "source": [ + "## Triangle and Quadrilateral Elements on a 2D Tire\n", + "\n", + "The tire is an annulus centred at the origin, with inner radius $R_i$ and outer radius $R_o$, adjusted to 2D from [[2]](#ref2)." + ] + }, + { + "cell_type": "markdown", + "id": "18-md-quartic-2d", + "metadata": {}, + "source": [ + "### Quartic Function\n", + "\n", + "$$\\hat{\\mathbf{u}} = A\\begin{pmatrix} x^4 + \\frac{2}{5}xy \\\\ \\frac{1}{10}y^4 - \\frac{2}{5}xy \\end{pmatrix}, \\qquad \\mathbf{s} = -\\frac{2A}{5}\\begin{pmatrix} 30(\\lambda + 2\\mu)x^2 - (\\lambda + \\mu) \\\\ 3(\\lambda + 2\\mu)y^2 + (\\lambda + \\mu) \\end{pmatrix}, \\qquad \\hat{\\boldsymbol{\\sigma}} = A\\begin{pmatrix} \\lambda\\theta + 2\\mu\\left(4x^3 + \\frac{2}{5}y\\right) & \\frac{2}{5}\\mu(x - y) \\\\ \\frac{2}{5}\\mu(x - y) & \\lambda\\theta + \\frac{4}{5}\\mu\\left(y^3 - x\\right) \\end{pmatrix}$$\n", + "\n", + "with $\\theta = 4x^3 + \\frac{2}{5}\\left(y^3 + y - x\\right)$.\n", + "\n", + "The boundary is the inner and outer rims of the annulus:\n", + "\n", + "- `inner` ($r = R_i$): Dirichlet, $\\mathbf{u} = \\hat{\\mathbf{u}}$\n", + "- `outer` ($r = R_o$): Neumann, $\\mathbf{t} = \\hat{\\boldsymbol{\\sigma}} \\cdot \\mathbf{n}$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19-code-quartic-2d", + "metadata": {}, + "outputs": [], + "source": [ + "# Case: 2D tire, quartic manufactured solution, triangles (left) and quadrilaterals (right)\n", + "tri = run(\"../cases/LinearSmallStrainFEMForceField/Quartic2D_Tri.json\")\n", + "quad = run(\"../cases/LinearSmallStrainFEMForceField/Quartic2D_Quad.json\")\n", + "figure, (triangles, quadrilaterals) = plt.subplots(1, 2, figsize=(14.4, 4.8))\n", + "tri.plot(triangles)\n", + "quad.plot(quadrilaterals)\n", + "\n", + "# Deformed shape at each level, coarse to fine\n", + "tri.plot_fields()\n", + "quad.plot_fields()" + ] + }, + { + "cell_type": "markdown", + "id": "20-md-3d", + "metadata": {}, + "source": [ + "## Tetrahedron and Hexahedron Elements on a 3D Beam\n", + "\n", + "The beam starts at the origin, with length $L$ along $x$, width $W$ along $y$ and height $H$ along $z$." + ] + }, + { + "cell_type": "markdown", + "id": "21-md-quadratic-3d", + "metadata": {}, + "source": [ + "### Quadratic Function\n", + "\n", + "$$\\hat{\\mathbf{u}} = A\\begin{pmatrix} x^2 + xy + xz \\\\ y^2 + yz + xy \\\\ z^2 + zx + yz \\end{pmatrix}, \\qquad \\mathbf{s} = -2A\\begin{pmatrix} 2\\lambda + 3\\mu \\\\ 2\\lambda + 3\\mu \\\\ 2\\lambda + 3\\mu \\end{pmatrix}, \\qquad \\hat{\\boldsymbol{\\sigma}} = A\\begin{pmatrix} 4\\lambda\\theta + 2\\mu(2x + y + z) & \\mu(x + y) & \\mu(x + z) \\\\ \\mu(x + y) & 4\\lambda\\theta + 2\\mu(x + 2y + z) & \\mu(y + z) \\\\ \\mu(x + z) & \\mu(y + z) & 4\\lambda\\theta + 2\\mu(x + y + 2z) \\end{pmatrix}$$\n", + "\n", + "with $\\theta = x + y + z$.\n", + "\n", + "- `left` ($x = 0$): $\\hat{u}_x = 0$, $\\hat{t}_y = -\\mu A y$, $\\hat{t}_z = -\\mu A z$\n", + "- `bottom` ($y = 0$): $\\hat{u}_y = 0$, $\\hat{t}_x = -\\mu A x$, $\\hat{t}_z = -\\mu A z$\n", + "- `front` ($z = 0$): $\\hat{u}_z = 0$, $\\hat{t}_x = -\\mu A x$, $\\hat{t}_y = -\\mu A y$\n", + "- `right` ($x = L$): $\\hat{\\mathbf{t}} = A\\,\\big(4\\lambda(L + y + z) + 2\\mu(2L + y + z),\\ \\mu(L + y),\\ \\mu(L + z)\\big)$\n", + "- `top` ($y = W$): $\\hat{\\mathbf{t}} = A\\,\\big(\\mu(x + W),\\ 4\\lambda(x + W + z) + 2\\mu(x + 2W + z),\\ \\mu(W + z)\\big)$\n", + "- `back` ($z = H$): $\\hat{\\mathbf{t}} = A\\,\\big(\\mu(x + H),\\ \\mu(y + H),\\ 4\\lambda(x + y + H) + 2\\mu(x + y + 2H)\\big)$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "22-code-quadratic-3d", + "metadata": {}, + "outputs": [], + "source": [ + "# Case: 3D box, quadratic manufactured solution, tetrahedra (left) and hexahedra (right)\n", + "tet = run(\"../cases/LinearSmallStrainFEMForceField/Quadratic3D_Tetra.json\")\n", + "hexa = run(\"../cases/LinearSmallStrainFEMForceField/Quadratic3D_Hexa.json\")\n", + "figure, (tetrahedra, hexahedra) = plt.subplots(1, 2, figsize=(14.4, 4.8))\n", + "tet.plot(tetrahedra)\n", + "hexa.plot(hexahedra)\n", + "\n", + "# Deformed shape at each level, coarse to fine\n", + "tet.plot_fields()\n", + "hexa.plot_fields()" + ] + }, + { + "cell_type": "markdown", + "id": "23-md-trigonometric-3d", + "metadata": {}, + "source": [ + "### Trigonometric Function\n", + "\n", + "$$\\hat{\\mathbf{u}} = A\\begin{pmatrix} \\sin(kx)\\cos(ky)\\cos(kz) \\\\ \\cos(kx)\\sin(ky)\\cos(kz) \\\\ \\cos(kx)\\cos(ky)\\sin(kz) \\end{pmatrix}, \\qquad \\mathbf{s} = 3k^2(\\lambda + 2\\mu)\\,\\hat{\\mathbf{u}}, \\qquad \\hat{\\boldsymbol{\\sigma}} = Ak\\begin{pmatrix} (3\\lambda + 2\\mu)\\,\\cos(kx)\\cos(ky)\\cos(kz) & -2\\mu\\,\\sin(kx)\\sin(ky)\\cos(kz) & -2\\mu\\,\\sin(kx)\\cos(ky)\\sin(kz) \\\\ -2\\mu\\,\\sin(kx)\\sin(ky)\\cos(kz) & (3\\lambda + 2\\mu)\\,\\cos(kx)\\cos(ky)\\cos(kz) & -2\\mu\\,\\cos(kx)\\sin(ky)\\sin(kz) \\\\ -2\\mu\\,\\sin(kx)\\cos(ky)\\sin(kz) & -2\\mu\\,\\cos(kx)\\sin(ky)\\sin(kz) & (3\\lambda + 2\\mu)\\,\\cos(kx)\\cos(ky)\\cos(kz) \\end{pmatrix}$$\n", + "\n", + "with $k = \\frac{2\\pi}{L}$, the same along $x$, $y$ and $z$ since the cases use a cube, $L = W = H$.\n", + "\n", + "- `left` ($x = 0$): $\\hat{u}_x = 0$, $\\hat{t}_y = \\hat{t}_z = 0$\n", + "- `bottom` ($y = 0$): $\\hat{u}_y = 0$, $\\hat{t}_x = \\hat{t}_z = 0$\n", + "- `front` ($z = 0$): $\\hat{u}_z = 0$, $\\hat{t}_x = \\hat{t}_y = 0$\n", + "- `right` ($x = L$): $\\hat{\\mathbf{t}} = Ak(3\\lambda + 2\\mu)\\,\\big(\\cos(ky)\\cos(kz),\\ 0,\\ 0\\big)$\n", + "- `top` ($y = L$): $\\hat{\\mathbf{t}} = Ak(3\\lambda + 2\\mu)\\,\\big(0,\\ \\cos(kx)\\cos(kz),\\ 0\\big)$\n", + "- `back` ($z = L$): $\\hat{\\mathbf{t}} = Ak(3\\lambda + 2\\mu)\\,\\big(0,\\ 0,\\ \\cos(kx)\\cos(ky)\\big)$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24-code-trigonometric-3d", + "metadata": {}, + "outputs": [], + "source": [ + "# Case: 3D box, trigonometric manufactured solution, tetrahedra (left) and hexahedra (right)\n", + "tet = run(\"../cases/LinearSmallStrainFEMForceField/Trigonometric3D_Tetra.json\")\n", + "hexa = run(\"../cases/LinearSmallStrainFEMForceField/Trigonometric3D_Hexa.json\")\n", + "figure, (tetrahedra, hexahedra) = plt.subplots(1, 2, figsize=(14.4, 4.8))\n", + "tet.plot(tetrahedra)\n", + "hexa.plot(hexahedra)\n", + "\n", + "# Deformed shape at each level, coarse to fine\n", + "tet.plot_fields()\n", + "hexa.plot_fields()" + ] + }, + { + "cell_type": "markdown", + "id": "25-md-3d-tire", + "metadata": {}, + "source": [ + "## Tetrahedron and Hexahedron Elements on a 3D Tire\n", + "\n", + "The tire is the annulus of the 2D tire, extruded along $z$ from $0$ to $H$." + ] + }, + { + "cell_type": "markdown", + "id": "26-md-quartic-3d", + "metadata": {}, + "source": [ + "### Quartic Function\n", + "\n", + "The manufactured solution is the one of [[2]](#ref2):\n", + "\n", + "$$\\hat{\\mathbf{u}} = A\\begin{pmatrix} x^4 - \\frac{2}{5}yz \\\\ -\\frac{1}{10}y^4 - \\frac{2}{5}xyz \\\\ z^4 - \\frac{2}{5}xy \\end{pmatrix}, \\qquad \\mathbf{s} = A\\begin{pmatrix} -12(\\lambda + 2\\mu)x^2 + \\frac{2}{5}(\\lambda + \\mu)z \\\\ \\frac{6}{5}(\\lambda + 2\\mu)y^2 \\\\ -12(\\lambda + 2\\mu)z^2 + \\frac{2}{5}(\\lambda + \\mu)x \\end{pmatrix}, \\qquad \\hat{\\boldsymbol{\\sigma}} = A\\begin{pmatrix} \\lambda\\theta + 8\\mu x^3 & -\\frac{2}{5}\\mu z(1 + y) & -\\frac{4}{5}\\mu y \\\\ -\\frac{2}{5}\\mu z(1 + y) & \\lambda\\theta - \\frac{4}{5}\\mu\\left(y^3 + xz\\right) & -\\frac{2}{5}\\mu x(1 + y) \\\\ -\\frac{4}{5}\\mu y & -\\frac{2}{5}\\mu x(1 + y) & \\lambda\\theta + 8\\mu z^3 \\end{pmatrix}$$\n", + "\n", + "with $\\theta = 4x^3 + 4z^3 - \\frac{2}{5}y^3 - \\frac{2}{5}xz$.\n", + "\n", + "The boundary is the inner and outer curved surfaces and the two flat faces:\n", + "\n", + "- `inner` ($r = R_i$): Dirichlet, $\\mathbf{u} = \\hat{\\mathbf{u}}$\n", + "- `outer` ($r = R_o$): Dirichlet, $\\mathbf{u} = \\hat{\\mathbf{u}}$\n", + "- `front` ($z = 0$): Neumann, $\\mathbf{t} = \\hat{\\boldsymbol{\\sigma}} \\cdot \\mathbf{n}$\n", + "- `back` ($z = H$): Neumann, $\\mathbf{t} = \\hat{\\boldsymbol{\\sigma}} \\cdot \\mathbf{n}$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "27-code-quartic-3d", + "metadata": {}, + "outputs": [], + "source": [ + "# Case: 3D tire, quartic manufactured solution, tetrahedra (left) and hexahedra (right)\n", + "tet = run(\"../cases/LinearSmallStrainFEMForceField/Quartic3D_Tetra.json\")\n", + "hexa = run(\"../cases/LinearSmallStrainFEMForceField/Quartic3D_Hexa.json\")\n", + "figure, (tetrahedra, hexahedra) = plt.subplots(1, 2, figsize=(14.4, 4.8))\n", + "tet.plot(tetrahedra)\n", + "hexa.plot(hexahedra)\n", + "\n", + "# Deformed shape at each level, coarse to fine\n", + "tet.plot_fields()\n", + "hexa.plot_fields()" + ] + }, + { + "cell_type": "markdown", + "id": "28-md-references", + "metadata": {}, + "source": [ + "\n", + "## References\n", + "\n", + "\n", + "**[1]** T. J. R. Hughes, *The Finite Element Method: Linear Static and Dynamic Finite Element Analysis*, Dover, Mineola, NY, 2000. ISBN 978-0-486-41181-1. — *Chapter 2 (the weak formulation of linear elasticity).*\n", + "\n", + "\n", + "**[2]** É. Chamberland, A. Fortin and M. Fortin, \"Comparison of the performance of some finite element discretizations for large deformation elasticity problems,\" *Computers and Structures*, 88(11–12):664–673, 2010." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/verification/manufactured/__init__.py b/verification/manufactured/__init__.py new file mode 100644 index 0000000..44adc87 --- /dev/null +++ b/verification/manufactured/__init__.py @@ -0,0 +1,5 @@ +from .base import ManufacturedSolution, TrigonometricSolution +from .problem import ManufacturedProblem +from .grid import (Quadratic1D, Quadratic2D, Quadratic3D, + Trigonometric1D, Trigonometric2D, Trigonometric3D) +from .tire import Quartic2D, Quartic3D diff --git a/verification/manufactured/base.py b/verification/manufactured/base.py new file mode 100644 index 0000000..4ec0f83 --- /dev/null +++ b/verification/manufactured/base.py @@ -0,0 +1,62 @@ +"""Manufactured solutions: exact displacements declared symbolically, and the regions where +their BCs apply. +""" + +from abc import ABC, abstractmethod + +import sympy as sp + + +class ManufacturedSolution(ABC): + """An exact displacement, declared symbolically, and the regions where its BCs apply. + + Stated without reference to a material; `ManufacturedProblem` pairs the two. + """ + + geometry: type # the geometry family this solution is stated on + dim: int # the PDE's own dimension, not the space it is solved in + + def __init__(self, config, geometry, traction_on=None): + if not isinstance(geometry, self.geometry): + raise ValueError(f"{type(self).__name__} is stated on a {self.geometry.__name__}, " + f"got a {type(geometry).__name__}") + if geometry.dim != self.dim: + raise ValueError(f"{type(self).__name__} is {self.dim}D, " + f"got a {geometry.dim}D {type(geometry).__name__}") + + self.config = config # the case's "solution" block; a subclass may read more from it + self.amplitude = config["amplitude"] # too large relative to the geometry inverts elements + self.parameters = geometry.named_parameters + self.regions = list(geometry.region_names) + self.traction_on = traction_on or {} # {region: directions} the case hands to traction + + @property + def prescribe_displacement_on(self): + """{region: fixed_directions}: every direction of every region, less the ones the case + handed to the traction. A region left with nothing fixed is dropped. + """ + prescribed = {} + for region in self.regions: + free = self.traction_on.get(region, [0] * self.dim) + fixed_directions = [1 - direction for direction in free] + if any(fixed_directions): + prescribed[region] = fixed_directions + return prescribed + + @abstractmethod + def displacement(self, coordinates): + """Exact displacement as sympy expressions, one per component of this solution's own dimension.""" + + +class TrigonometricSolution(ManufacturedSolution): + """A solution built from sines and cosines fitted to the geometry's extent.""" + + def __init__(self, config, geometry, traction_on=None): + super().__init__(config, geometry, traction_on) + self.periods = config.get("periods", 1) # more periods need a finer mesh to resolve + + def wavenumber(self, parameter): + """The wavenumber fitting `self.periods` full periods across one of the geometry's + characteristic parameters, e.g. 2 pi / length for a single period. + """ + return 2 * sp.pi * sp.Rational(str(self.periods)) / sp.Rational(str(self.parameters[parameter])) diff --git a/verification/manufactured/grid.py b/verification/manufactured/grid.py new file mode 100644 index 0000000..b43d346 --- /dev/null +++ b/verification/manufactured/grid.py @@ -0,0 +1,97 @@ +"""Manufactured solutions stated on the grid geometry.""" + +import sympy as sp + +from common.geometry import Grid + +from .base import ManufacturedSolution, TrigonometricSolution + + +class Quadratic1D(ManufacturedSolution): + """u(x) = [A x^2]: body force constant, so the nodal source is the source; clamped at 'left'.""" + + geometry = Grid + dim = 1 + + def displacement(self, coordinates): + return [self.amplitude * coordinates[0] ** 2] + + +class Quadratic2D(ManufacturedSolution): + """u(x, y) = [A (x^2 + x y), A (y^2 + x y)] + + Body force constant and traction linear, so both are nodally exact. + """ + + geometry = Grid + dim = 2 + + def displacement(self, coordinates): + x, y = coordinates[0], coordinates[1] + return [self.amplitude * (x**2 + x * y), self.amplitude * (y**2 + x * y)] + + +class Quadratic3D(ManufacturedSolution): + """u(x, y, z) = [A (x^2 + x y + x z), A (y^2 + y z + x y), A (z^2 + z x + y z)] + + Body force constant and traction linear, so both are nodally exact. + """ + + geometry = Grid + dim = 3 + + def displacement(self, coordinates): + x, y, z = coordinates + return [self.amplitude * (x**2 + x * y + x * z), + self.amplitude * (y**2 + y * z + x * y), + self.amplitude * (z**2 + z * x + y * z)] + + +class Trigonometric1D(TrigonometricSolution): + """u(x) = [A sin(k x)], k = 2 pi / L: prescribed (clamped) at 'left', traction at 'right'.""" + + geometry = Grid + dim = 1 + + def displacement(self, coordinates): + return [self.amplitude * sp.sin(self.wavenumber("length") * coordinates[0])] + + +class Trigonometric2D(TrigonometricSolution): + """u(x, y) = [A sin(kx x) cos(ky y), + A cos(kx x) sin(ky y)] + kx = 2 pi / L, ky = 2 pi / W + + u_x fixed on the left face, u_y on the bottom face, traction elsewhere. + """ + + geometry = Grid + dim = 2 + + def displacement(self, coordinates): + kx, ky = self.wavenumber("length"), self.wavenumber("width") + x, y = coordinates[0], coordinates[1] + return [self.amplitude * sp.sin(kx * x) * sp.cos(ky * y), + self.amplitude * sp.cos(kx * x) * sp.sin(ky * y)] + + +class Trigonometric3D(TrigonometricSolution): + """u(x, y, z) = [A sin(kx x) cos(ky y) cos(kz z), + A cos(kx x) sin(ky y) cos(kz z), + A cos(kx x) cos(ky y) sin(kz z)] + kx = 2 pi / L, ky = 2 pi / W, kz = 2 pi / H + + Each component fixed on the low face of its own axis, traction elsewhere. + """ + + geometry = Grid + dim = 3 + + def displacement(self, coordinates): + kx, ky, kz = (self.wavenumber("length"), self.wavenumber("width"), + self.wavenumber("height")) + x, y, z = coordinates + return [self.amplitude * sp.sin(kx * x) * sp.cos(ky * y) * sp.cos(kz * z), + self.amplitude * sp.cos(kx * x) * sp.sin(ky * y) * sp.cos(kz * z), + self.amplitude * sp.cos(kx * x) * sp.cos(ky * y) * sp.sin(kz * z)] + diff --git a/verification/manufactured/materials.py b/verification/manufactured/materials.py new file mode 100644 index 0000000..d804836 --- /dev/null +++ b/verification/manufactured/materials.py @@ -0,0 +1,131 @@ +"""Constitutive laws: a material declares psi(gradient); its stress and tangent are differentiated from that.""" + +from abc import ABC, abstractmethod + +import numpy as np +import sympy as sp + + +def _toLameParameters1D(youngModulus, poissonRatio): + """(mu, lambda) at d = 1, where lambda + 2 mu = E collapses Hooke to sigma = E eps.""" + return 0.5 * youngModulus, 0.0 + + +def _toLameParameters2D(youngModulus, poissonRatio): + """(mu, lambda) at d = 2, in plane stress: the convention SOFA's force fields assume.""" + return (youngModulus / (2 * (1 + poissonRatio)), + youngModulus * poissonRatio / (1 - poissonRatio ** 2)) + + +def _toLameParameters3D(youngModulus, poissonRatio): + """(mu, lambda) at d = 3.""" + return (youngModulus / (2 * (1 + poissonRatio)), + youngModulus * poissonRatio / ((1 + poissonRatio) * (1 - 2 * poissonRatio))) + + +# dim -> Young/Poisson -> (mu, lambda) +_toLame = {1: _toLameParameters1D, + 2: _toLameParameters2D, + 3: _toLameParameters3D} + + +def _tensor(dimensions, label): + """A d x d matrix of independent symbols, named "