From 94d0ae80608d42474c2f9b2a2d0a167caa901bc6 Mon Sep 17 00:00:00 2001 From: srosenbu Date: Fri, 19 Jan 2024 13:01:39 +0100 Subject: [PATCH 01/19] Fix type hints --- src/fenics_constitutive/interfaces.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/fenics_constitutive/interfaces.py b/src/fenics_constitutive/interfaces.py index 58f0a37..1f2aaa1 100644 --- a/src/fenics_constitutive/interfaces.py +++ b/src/fenics_constitutive/interfaces.py @@ -53,7 +53,7 @@ def evaluate( grad_del_u: np.ndarray, mandel_stress: np.ndarray, tangent: np.ndarray, - history: np.ndarray | dict[str, np.ndarray], + history: np.ndarray | dict[str, np.ndarray] | None, ) -> None: """Evaluate the constitutive model and overwrite the stress, tangent and history. @@ -107,7 +107,7 @@ def geometric_dim(self) -> int: return self.constraint.geometric_dim() @abstractproperty - def history_dim(self) -> int | dict[str, int | tuple[int, int]]: + def history_dim(self) -> int | dict[str, int | tuple[int, int]] | None: """ The dimensions of history variable(s). This is needed to tell the solver which quadrature spaces or arrays to build. If all history variables are stored in a single From ec5fb68a9c6311d82f1ef06cace87cc4b2053a04 Mon Sep 17 00:00:00 2001 From: srosenbu Date: Mon, 22 Jan 2024 13:28:50 +0100 Subject: [PATCH 02/19] linear elasticity for 3d --- .clang-format | 26 ------- .../linear_elasticity_model.py | 53 +++++++++++++ src/fenics_constitutive/__init__.py | 1 + src/fenics_constitutive/_version.py | 4 +- src/fenics_constitutive/conversion.py | 77 +++++++++++++++++++ tests/test_conversions.py | 35 +++++++++ 6 files changed, 168 insertions(+), 28 deletions(-) delete mode 100644 .clang-format create mode 100644 examples/linear_elasticity/linear_elasticity_model.py create mode 100644 src/fenics_constitutive/conversion.py create mode 100644 tests/test_conversions.py diff --git a/.clang-format b/.clang-format deleted file mode 100644 index 4e3cd54..0000000 --- a/.clang-format +++ /dev/null @@ -1,26 +0,0 @@ -BasedOnStyle: LLVM -Language: Cpp - -AlignAfterOpenBracket: Align -AlignOperands: true -AlignConsecutiveAssignments: false -AlignConsecutiveDeclarations: false -AlignTrailingComments: false -AlwaysBreakTemplateDeclarations: true -AllowShortFunctionsOnASingleLine: true -AccessModifierOffset: -4 -BreakBeforeBraces: Allman -ColumnLimit: 120 -AllowShortIfStatementsOnASingleLine: false -AllowShortBlocksOnASingleLine: true -AllowShortFunctionsOnASingleLine: false -ContinuationIndentWidth: 8 -DerivePointerAlignment: false -IndentWidth: 4 -PointerAlignment: Left -ReflowComments: true -MaxEmptyLinesToKeep: 2 -BreakConstructorInitializersBeforeComma: true -SpaceBeforeParens: ControlStatements -SortIncludes: false -Standard: Cpp11 diff --git a/examples/linear_elasticity/linear_elasticity_model.py b/examples/linear_elasticity/linear_elasticity_model.py new file mode 100644 index 0000000..7d678ab --- /dev/null +++ b/examples/linear_elasticity/linear_elasticity_model.py @@ -0,0 +1,53 @@ +from fenics_constitutive import Constraint, IncrSmallStrainModel, strain_from_grad_u +import numpy as np + + +class LinearElasticityModel(IncrSmallStrainModel): + def __init__(self, E: float, nu: float, constraint: Constraint): + self._constraint = constraint + mu = E / (2 * (1 + nu)) + lam = E * nu / ((1 + nu) * (1 - 2 * nu)) + match constraint: + case Constraint.FULL: + self.D = np.array( + [ + [2.0 * mu + lam, lam, lam, 0.0, 0.0, 0.0], + [lam, 2.0 * mu + lam, lam, 0.0, 0.0, 0.0], + [lam, lam, 2.0 * mu + lam, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, mu, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, mu, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, mu], + ] + ) + case _: + raise NotImplementedError("Only full constraint implemented") + + def evaluate( + self, + del_t: float, + grad_del_u: np.ndarray, + mandel_stress: np.ndarray, + tangent: np.ndarray, + history: np.ndarray | dict[str, np.ndarray] | None, + ) -> None: + assert ( + grad_del_u.size / (self.geometric_dim**2) + == mandel_stress.size / self.stress_strain_dim + == tangent.size / (self.stress_strain_dim**2) + ) + n_gauss = grad_del_u.shape / (self.geometric_dim**2) + + strain_increment = strain_from_grad_u(grad_del_u, self.constraint) + mandel_stress += strain_increment.reshape(-1, self.stress_strain_dim) @ self.D + tangent[:] = np.tile(self.D, n_gauss) + + @property + def constraint(self) -> Constraint: + return self._constraint + + @property + def history_dim(self) -> None: + return None + + def update(self) -> None: + pass \ No newline at end of file diff --git a/src/fenics_constitutive/__init__.py b/src/fenics_constitutive/__init__.py index 562dd05..a337942 100644 --- a/src/fenics_constitutive/__init__.py +++ b/src/fenics_constitutive/__init__.py @@ -8,6 +8,7 @@ from __future__ import annotations from .interfaces import * +from .conversion import * from ._version import version as __version__ __all__ = ["__version__"] diff --git a/src/fenics_constitutive/_version.py b/src/fenics_constitutive/_version.py index 98515de..5c2bafa 100644 --- a/src/fenics_constitutive/_version.py +++ b/src/fenics_constitutive/_version.py @@ -12,5 +12,5 @@ __version_tuple__: VERSION_TUPLE version_tuple: VERSION_TUPLE -__version__ = version = '0.1.dev104+g7b69162' -__version_tuple__ = version_tuple = (0, 1, 'dev104', 'g7b69162') +__version__ = version = '0.1.dev110+ge118955' +__version_tuple__ = version_tuple = (0, 1, 'dev110', 'ge118955') diff --git a/src/fenics_constitutive/conversion.py b/src/fenics_constitutive/conversion.py new file mode 100644 index 0000000..6161917 --- /dev/null +++ b/src/fenics_constitutive/conversion.py @@ -0,0 +1,77 @@ +import numpy as np +from .interfaces import Constraint + +__all__ = ["strain_from_grad_u"] + +def strain_from_grad_u(grad_u: np.ndarray, constraint: Constraint) -> np.ndarray: + """Compute the Mandel-strain from the gradient of displacement. + + Parameters: + grad_u: Gradient of displacement field. + + Returns: + Strain tensor. + """ + strain_dim = constraint.stress_strain_dim() + n_gauss = int(grad_u.size / (constraint.geometric_dim()**2)) + strain = np.zeros(strain_dim*n_gauss) + grad_u_view = grad_u.reshape(-1, constraint.geometric_dim()**2) + strain_view = strain.reshape(-1, strain_dim) + + match constraint: + case Constraint.UNIAXIAL_STRAIN: + strain_view[:, 0] = grad_u_view[:, 0] + case Constraint.UNIAXIAL_STRESS: + strain_view[:,0] = grad_u_view[:, 0] + case Constraint.PLANE_STRAIN: + """ + Full tensor: + + 0 1 + 2 3 + + Mandel vector: + f = 1 / sqrt(2) + 0 3 "0" f*(1+2) + """ + strain_view[:,0] = grad_u_view[:,0] + strain_view[:,1] = grad_u_view[:,3] + strain_view[:,2] = 0.0 + strain_view[:,3] = 1/2**0.5 * (grad_u_view[:,1] + grad_u_view[:,2]) + case Constraint.PLANE_STRESS: + """ + Full tensor: + + 0 1 + 2 3 + + Mandel vector: + f = 1 / sqrt(2) + 0 3 "0" f*(1+2) + """ + strain_view[:,0] = grad_u_view[:,0] + strain_view[:,1] = grad_u_view[:,3] + strain_view[:,2] = 0.0 + strain_view[:,3] = 1/2**0.5 * (grad_u_view[:,1] + grad_u_view[:,2]) + case Constraint.FULL: + """ + Full tensor: + + 0 1 2 + 3 4 5 + 6 7 8 + + Mandel vector: + f = 1 / sqrt(2) + 0 4 8 f*(1+3) f*(5+7) f*(2+6) + """ + strain_view[:,0] = grad_u_view[:,0] + strain_view[:,1] = grad_u_view[:,4] + strain_view[:,2] = grad_u_view[:,8] + strain_view[:,3] = 1/2**0.5 * (grad_u_view[:,1] + grad_u_view[:,3]) + strain_view[:,4] = 1/2**0.5 * (grad_u_view[:,5] + grad_u_view[:,7]) + strain_view[:,5] = 1/2**0.5 * (grad_u_view[:,2] + grad_u_view[:,6]) + case _: + raise NotImplementedError("Only full constraint implemented") + return strain + \ No newline at end of file diff --git a/tests/test_conversions.py b/tests/test_conversions.py new file mode 100644 index 0000000..aaf3af4 --- /dev/null +++ b/tests/test_conversions.py @@ -0,0 +1,35 @@ +from fenics_constitutive import Constraint, strain_from_grad_u +import numpy as np + + +def test_strain_from_grad_u(): + grad_u = np.array([[1.0]]) + constraint = Constraint.UNIAXIAL_STRAIN + strain = strain_from_grad_u(grad_u, constraint) + assert np.allclose(strain, np.array([1.0])) + constraint = Constraint.UNIAXIAL_STRESS + strain = strain_from_grad_u(grad_u, constraint) + assert np.allclose(strain, np.array([1.0])) + grad_u = np.array([[1.0, 2.0], [3.0, 4.0]]) + constraint = Constraint.PLANE_STRAIN + strain = strain_from_grad_u(grad_u, constraint) + assert np.allclose(strain, np.array([1.0, 4.0, 0.0, 0.5 * (4.0 + 1.0) * 2**0.5])) + constraint = Constraint.PLANE_STRESS + strain = strain_from_grad_u(grad_u, constraint) + assert np.allclose(strain, np.array([1.0, 4.0, 0.0, 0.5 * (4.0 + 1.0) * 2**0.5])) + grad_u = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]) + constraint = Constraint.FULL + strain = strain_from_grad_u(grad_u, constraint) + assert np.allclose( + strain, + np.array( + [ + 1.0, + 5.0, + 9.0, + 0.5 * (2.0 + 4.0) * 2**0.5, + 0.5 * (6.0 + 8.0) * 2**0.5, + 0.5 * (3.0 + 7.0) * 2**0.5, + ] + ), + ) From e5fd7d11ab429aa8d830cc3108068a3bff6aa43c Mon Sep 17 00:00:00 2001 From: srosenbu Date: Wed, 31 Jan 2024 16:45:03 +0100 Subject: [PATCH 03/19] All constraints --- .../linear_elasticity_model.py | 51 ++++++++++++++++--- 1 file changed, 43 insertions(+), 8 deletions(-) diff --git a/examples/linear_elasticity/linear_elasticity_model.py b/examples/linear_elasticity/linear_elasticity_model.py index 7d678ab..fdd3a79 100644 --- a/examples/linear_elasticity/linear_elasticity_model.py +++ b/examples/linear_elasticity/linear_elasticity_model.py @@ -9,18 +9,53 @@ def __init__(self, E: float, nu: float, constraint: Constraint): lam = E * nu / ((1 + nu) * (1 - 2 * nu)) match constraint: case Constraint.FULL: + # see https://en.wikipedia.org/wiki/Hooke%27s_law self.D = np.array( [ [2.0 * mu + lam, lam, lam, 0.0, 0.0, 0.0], [lam, 2.0 * mu + lam, lam, 0.0, 0.0, 0.0], [lam, lam, 2.0 * mu + lam, 0.0, 0.0, 0.0], - [0.0, 0.0, 0.0, mu, 0.0, 0.0], - [0.0, 0.0, 0.0, 0.0, mu, 0.0], - [0.0, 0.0, 0.0, 0.0, 0.0, mu], + [0.0, 0.0, 0.0, 2.0 * mu, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 2.0 * mu, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 2.0 * mu], ] ) + case Constraint.PLANE_STRAIN: + # We assert that the strain is being provided with 0 in the z-direction + # see https://en.wikipedia.org/wiki/Hooke%27s_law + self.D = np.array( + [ + [2.0 * mu + lam, lam, 0.0, 0.0], + [lam, 2.0 * mu + lam, 0.0, 0.0], + [lam, lam, 2.0 * mu + lam, 0.0], + [0.0, 0.0, 0.0, 2.0 * mu], + ] + ) + case Constraint.PLANE_STRESS: + # We assert that the strain is being provided with 0 in the z-direction + # This matrix just multiplies the z component by 1.0 + # see https://en.wikipedia.org/wiki/Hooke%27s_law + self.D = ( + E + / (1 - nu**2.0) + * np.array( + [ + [1.0, nu, 0.0, 0.0], + [nu, 1.0, 0.0, 0.0], + [0.0, 0.0, 1.0, 0.0], + [0.0, 0.0, 0.0, (1.0 - nu)], + ] + ) + ) + case Constraint.UNIAXIAL_STRAIN: + # see https://csmbrannon.net/2012/08/02/distinction-between-uniaxial-stress-and-uniaxial-strain/ + C = E * (1.0 - nu) / ((1.0 + nu) * (1.0 - 2.0 * nu)) + self.D = np.array([[C]]) + case Constraint.UNIAXIAL_STRESS: + # see https://csmbrannon.net/2012/08/02/distinction-between-uniaxial-stress-and-uniaxial-strain/ + self.D = np.array([[E]]) case _: - raise NotImplementedError("Only full constraint implemented") + raise NotImplementedError("Constraint not implemented") def evaluate( self, @@ -40,14 +75,14 @@ def evaluate( strain_increment = strain_from_grad_u(grad_del_u, self.constraint) mandel_stress += strain_increment.reshape(-1, self.stress_strain_dim) @ self.D tangent[:] = np.tile(self.D, n_gauss) - + @property def constraint(self) -> Constraint: return self._constraint - + @property def history_dim(self) -> None: return None - + def update(self) -> None: - pass \ No newline at end of file + pass From b567c9bc5a35e86f704f59af82a1d4d9518b8ec1 Mon Sep 17 00:00:00 2001 From: srosenbu Date: Wed, 31 Jan 2024 16:48:56 +0100 Subject: [PATCH 04/19] parameters as dict --- examples/linear_elasticity/linear_elasticity_model.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/examples/linear_elasticity/linear_elasticity_model.py b/examples/linear_elasticity/linear_elasticity_model.py index fdd3a79..761239f 100644 --- a/examples/linear_elasticity/linear_elasticity_model.py +++ b/examples/linear_elasticity/linear_elasticity_model.py @@ -3,10 +3,12 @@ class LinearElasticityModel(IncrSmallStrainModel): - def __init__(self, E: float, nu: float, constraint: Constraint): + def __init__(self, parameters: dict[str, float], constraint: Constraint): self._constraint = constraint - mu = E / (2 * (1 + nu)) - lam = E * nu / ((1 + nu) * (1 - 2 * nu)) + E = parameters["E"] + nu = parameters["nu"] + mu = E / (2.0 * (1.0 + nu)) + lam = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu)) match constraint: case Constraint.FULL: # see https://en.wikipedia.org/wiki/Hooke%27s_law From a2a84c7926b4c4e402ba0ac4456c5acfcbe8e948 Mon Sep 17 00:00:00 2001 From: srosenbu Date: Tue, 19 Mar 2024 16:00:24 +0100 Subject: [PATCH 05/19] Work on solver and mappings --- .../linear_elasticity_model.py | 7 +- src/fenics_constitutive/__init__.py | 2 + src/fenics_constitutive/maps.py | 139 +++++++++++ src/fenics_constitutive/solver.py | 220 ++++++++++++++++++ tests/test_maps.py | 102 ++++++++ 5 files changed, 467 insertions(+), 3 deletions(-) create mode 100644 src/fenics_constitutive/maps.py create mode 100644 src/fenics_constitutive/solver.py create mode 100644 tests/test_maps.py diff --git a/examples/linear_elasticity/linear_elasticity_model.py b/examples/linear_elasticity/linear_elasticity_model.py index 761239f..4bbb00e 100644 --- a/examples/linear_elasticity/linear_elasticity_model.py +++ b/examples/linear_elasticity/linear_elasticity_model.py @@ -34,8 +34,9 @@ def __init__(self, parameters: dict[str, float], constraint: Constraint): ] ) case Constraint.PLANE_STRESS: - # We assert that the strain is being provided with 0 in the z-direction - # This matrix just multiplies the z component by 1.0 + # We do not make any assumptions about strain in the z-direction + # This matrix just multiplies the z component by 0.0 which results + # in a plane stress state # see https://en.wikipedia.org/wiki/Hooke%27s_law self.D = ( E @@ -44,7 +45,7 @@ def __init__(self, parameters: dict[str, float], constraint: Constraint): [ [1.0, nu, 0.0, 0.0], [nu, 1.0, 0.0, 0.0], - [0.0, 0.0, 1.0, 0.0], + [0.0, 0.0, 0.0, 0.0], [0.0, 0.0, 0.0, (1.0 - nu)], ] ) diff --git a/src/fenics_constitutive/__init__.py b/src/fenics_constitutive/__init__.py index a41c97c..f001bb6 100644 --- a/src/fenics_constitutive/__init__.py +++ b/src/fenics_constitutive/__init__.py @@ -10,5 +10,7 @@ from ._version import version as __version__ from .interfaces import * from .stress_strain import * +from .solver import * +from .maps import * __all__ = ["__version__"] diff --git a/src/fenics_constitutive/maps.py b/src/fenics_constitutive/maps.py new file mode 100644 index 0000000..748f098 --- /dev/null +++ b/src/fenics_constitutive/maps.py @@ -0,0 +1,139 @@ +import dolfinx as df +import numpy as np +from functools import reduce +from dataclasses import dataclass + +__all__ = ["SubSpaceMap", "_build_view", "build_subspace_map"] + + +@dataclass +class SubSpaceMap: + parent: np.ndarray + child: np.ndarray + sub_mesh: df.mesh.Mesh + parent_mesh: df.mesh.Mesh + + def map_to_parent(self, sub: df.fem.Function, parent: df.fem.Function) -> None: + assert sub.function_space.mesh == self.sub_mesh, "Subspace mesh does not match" + assert ( + parent.function_space.mesh == self.parent_mesh + ), "Parent mesh does not match" + assert sub.ufl_shape == parent.ufl_shape, "Shapes do not match" + # print(sub.ufl_element().value_size(), "sub.ufl_element().value_size()") + size = sub.ufl_element().value_size() # reduce(lambda x, y: x * y, sub.ufl_shape) if not isinstance(sub.ufl_shape, int) else sub.ufl_shape + + parent_array = parent.x.array.reshape(-1, size) + sub_array = sub.x.array.reshape(-1, size) + parent_array[self.parent] = sub_array[self.child] + + def map_to_child(self, parent: df.fem.Function, sub: df.fem.Function) -> None: + assert sub.function_space.mesh == self.sub_mesh, "Subspace mesh does not match" + assert ( + parent.function_space.mesh == self.parent_mesh + ), "Parent mesh does not match" + assert sub.ufl_shape == parent.ufl_shape, "Shapes do not match" + + # size = reduce(lambda x, y: x * y, sub.ufl_shape) if not isinstance(sub.ufl_shape, int) else sub.ufl_shape + size = sub.ufl_element().value_size() # reduce(lambda x, y: x * y, sub.ufl_shape) if not isinstance(sub.ufl_shape, int) else sub.ufl_shape + + parent_array = parent.x.array.reshape(-1, size) + sub_array = sub.x.array.reshape(-1, size) + parent_array[self.parent] = sub_array[self.child] + + +def build_subspace_map( + cells: np.ndarray, V: df.fem.FunctionSpace, return_subspace=False +) -> ( + tuple[SubSpaceMap, df.mesh.Mesh] + | tuple[SubSpaceMap, df.mesh.Mesh, df.fem.FunctionSpace] +): + mesh = V.mesh + submesh, cell_map, _, _ = df.mesh.create_submesh(mesh, mesh.topology.dim, cells) + fe = V.ufl_element() + V_sub = df.fem.FunctionSpace(submesh, fe) + + submesh = V_sub.mesh + view_parent = [] + view_child = [] + + num_sub_cells = submesh.topology.index_map(submesh.topology.dim).size_local + for cell in range(num_sub_cells): + view_child.append(V_sub.dofmap.cell_dofs(cell)) + view_parent.append(V.dofmap.cell_dofs(cell_map[cell])) + + if view_child: + map = SubSpaceMap( + parent=np.hstack(view_parent), + child=np.hstack(view_child), + sub_mesh=submesh, + parent_mesh=V.mesh, + ) + else: + map = SubSpaceMap( + parent=np.array([], dtype=np.int32), + child=np.array([], dtype=np.int32), + sub_mesh=submesh, + parent_mesh=V.mesh, + ) + if return_subspace: + return map, submesh, V_sub + else: + del V_sub + return map, submesh + + +def _build_view( + cells: np.ndarray, V: df.fem.FunctionSpace +) -> tuple[SubSpaceMap, df.fem.FunctionSpace]: + mesh = V.mesh + submesh, cell_map, _, _ = df.mesh.create_submesh(mesh, mesh.topology.dim, cells) + fe = V.ufl_element() + V_sub = df.fem.FunctionSpace(submesh, fe) + + submesh = V_sub.mesh + view_parent = [] + view_child = [] + + num_sub_cells = submesh.topology.index_map(submesh.topology.dim).size_local + for cell in range(num_sub_cells): + view_child.append(V_sub.dofmap.cell_dofs(cell)) + view_parent.append(V.dofmap.cell_dofs(cell_map[cell])) + if view_child: + return ( + SubSpaceMap( + parent=np.hstack(view_parent), + child=np.hstack(view_child), + sub_mesh=submesh, + parent_mesh=V.mesh, + ), + V_sub, + ) + else: + # it may be that a process does not own any of the cells in the submesh + return ( + SubSpaceMap( + parent=np.array([], dtype=np.int32), + child=np.array([], dtype=np.int32), + sub_mesh=submesh, + parent_mesh=V.mesh, + ), + V_sub, + ) + return ( + np.array([], dtype=np.int32), + np.array([], dtype=np.int32), + V_sub, + ) + # if view_child: + # return ( + # np.hstack(view_parent), + # np.hstack(view_child), + # V_sub, + # ) + # else: + # # it may be that a process does not own any of the cells in the submesh + # return ( + # np.array([], dtype=np.int32), + # np.array([], dtype=np.int32), + # V_sub, + # ) diff --git a/src/fenics_constitutive/solver.py b/src/fenics_constitutive/solver.py new file mode 100644 index 0000000..9ba3205 --- /dev/null +++ b/src/fenics_constitutive/solver.py @@ -0,0 +1,220 @@ +from dataclasses import dataclass +from typing import Any +import basix +import dolfinx as df +import numpy as np +import ufl +from petsc4py import PETSc +from .interfaces import IncrSmallStrainModel +from .stress_strain import ufl_mandel_strain + +from .maps import SubSpaceMap, build_subspace_map + +def build_history( + law: IncrSmallStrainModel, submesh: df.mesh.Mesh, q_degree: int +) -> df.fem.Function | dict[str, df.fem.Function] | None: + match law.history_dim: + case int(): + Qh = ufl.VectorElement( + "Quadrature", + submesh.ufl_cell(), + q_degree, + quad_scheme="default", + dim=law.history_dim, + ) + history_space = df.fem.FunctionSpace(submesh, Qh) + history = df.fem.Function(history_space) + case None: + history = None + case dict(): + history = {} + for key, value in law.history_dim.items(): + if isinstance(value, int): + Qh = ufl.VectorElement( + "Quadrature", + submesh.ufl_cell(), + q_degree, + quad_scheme="default", + dim=value, + ) + elif isinstance(value, tuple): + Qh = ufl.TensorElement( + "Quadrature", + submesh.ufl_cell(), + q_degree, + quad_scheme="default", + shape=value, + ) + history_space = df.fem.FunctionSpace(submesh, Qh) + history[key] = df.fem.Function(history_space) + return history + + + +class IncrMechanicsProblem(df.fem.petsc.NonlinearProblem): + def __init__( + self, + # laws = [(law_1, cells_1), (law_2, cells_2), ...] + laws: list[tuple[IncrSmallStrainModel, np.ndarray]], + # u0 is basically a history variable + u0: df.fem.Function, + bcs: list[df.fem.DirichletBCMetaClass], + q_degree: int = 1, + form_compiler_options: dict = {}, + jit_options: dict = {}, + ): + mesh = u0.function_space.mesh + map_c = mesh.topology.index_map(mesh.topology.dim) + num_cells = map_c.size_local + map_c.num_ghosts + self.cells = np.arange(0, num_cells, dtype=np.int32) + + self.gdim = mesh.ufl_cell().geometric_dimension() + + constraint = laws[0][0].constraint + assert all( + law[0].constraint == constraint for law in laws + ), "All laws must have the same constraint" + + QVe = ufl.VectorElement( + "Quadrature", + mesh.ufl_cell(), + q_degree, + quad_scheme="default", + dim=constraint.stress_strain_dim(), + ) + QTe = ufl.TensorElement( + "Quadrature", + mesh.ufl_cell(), + q_degree, + quad_scheme="default", + shape=(constraint.stress_strain_dim(), constraint.stress_strain_dim()), + ) + QV = df.fem.FunctionSpace(mesh, QVe) + QT = df.fem.FunctionSpace(mesh, QTe) + + self.laws = [] + # TODO: We can probably reduce to one map per submesh + self.submesh_maps = [] + #self.QV_views = [] + #self.Qh_views = [] + #self.QT_views = [] + + self._strain = [] + self._stress = [] + self._history_0 = [] + self._history_1 = [] + self._tangent = [] + + self._time = 0.0 # time at the end of the increment + + with df.common.Timer("submeshes-and-data-structures"): + for law, cells in laws: + self.laws.append((law, cells)) + + + # ### submesh and subspace for strain, stress + subspace_map, submesh, QV_subspace = build_subspace_map(cells, QV, return_subspace=True) + self.submesh_maps.append(subspace_map) + self._strain.append(df.fem.Function(QV_subspace)) + self._stress.append(df.fem.Function(QV_subspace)) + + # ### submesh and subspace for tanget + QT_subspace = df.fem.FunctionSpace(submesh, QTe) + self._tangent.append(df.fem.Function(QT_subspace)) + + # ### submesh and subspace for history + history_0, Qh_views = build_history(law, submesh, q_degree) + history_1 = ( + {key: fn.copy() for key, fn in history_0.items()} + if isinstance(history_0, dict) + else history_0 + ) + self._history_0.append(history_0) + self._history_1.append(history_1) + + self.stress_0 = df.fem.Function(QV) + self.stress_1 = df.fem.Function(QV) + self.tangent = df.fem.Function(QT) + + u_, du = ufl.TestFunction(u0.function_space), ufl.TrialFunction(u0.function_space) + + self.metadata = {"quadrature_degree": q_degree, "quadrature_scheme": "default"} + self.dxm = ufl.dx(metadata=self.metadata) + + self.R_form = ufl.inner(ufl_mandel_strain(u_, constraint), self.stress) * self.dxm + self.dR_form = ufl.inner(ufl_mandel_strain(du, constraint), ufl.dot(self.tangent, ufl_mandel_strain(u_, constraint))) * self.dxm + + self._u0 = u0 + self._bcs = bcs + self._form_compiler_options = form_compiler_options + self._jit_options = jit_options + + basix_celltype = getattr(basix.CellType, mesh.topology.cell_type.name) + self.q_points, _ = basix.make_quadrature(basix_celltype, q_degree) + + # ### Note on JIT compilation of UFL forms + # if super().__init__(R, u, bcs, dR) is called within ElasticityProblem.__init__ + # the user cannot add Neumann BCs. + # Therefore, the compilation (i.e. call to super().__init__()) is done when + # df.nls.petsc.NewtonSolver is initialized. + # df.nls.petsc.NewtonSolver will call + # self._A = fem.petsc.create_matrix(problem.a) and hence it would suffice + # to override the property ``a`` of NonlinearProblem. + + @property + def a(self) -> df.fem.FormMetaClass: + """Compiled bilinear form (the Jacobian form)""" + + if not hasattr(self, "_a"): + # ensure compilation of UFL forms + super().__init__( + self.R_form, + self._u, + self._bcs, + self.dR_form, + form_compiler_options=self._form_compiler_options, + jit_options=self._jit_options, + ) + + return self._a + + def form(self, x: PETSc.Vec): + """This function is called before the residual or Jacobian is + computed. This is usually used to update ghost values. + Parameters + ---------- + x + The vector containing the latest solution + """ + super().form(x) + #TODO: look if the creation of the expression slows things down + strain_expr = df.fem.Expression(ufl_mandel_strain(x - self._u0, self.constraint), self.q_points) + + for k, (law, cells) in enumerate(self.laws): + with df.common.Timer("strain_evaluation"): + strain_expr.eval( + cells, self._strain[k].x.array.reshape(cells.size, -1) + ) + + with df.common.Timer("stress_evaluation"): + self.submesh_maps[k].map_to_child(self.stress_0, self._stress[k]) + self._history_1[k].x.array[:] = self._history_0[k].x.array + law.evaluate( + self._time, + self._strain[k].x.array, + self._stress[k].x.array, + self._tangent[k].x.array, + self._history_1[k].x.array, # history, + ) + + + with df.common.Timer("stress-local-to-global"): + self.submesh_maps[k].map_to_parent(self._stress[k], self.stress_1) + self.submesh_maps[k].map_to_parent(self._tangent[k], self.tangent) + + self.stress.x.scatter_forward() + self.tangent.x.scatter_forward() + + def update(self): + pass + diff --git a/tests/test_maps.py b/tests/test_maps.py new file mode 100644 index 0000000..836bc76 --- /dev/null +++ b/tests/test_maps.py @@ -0,0 +1,102 @@ +import dolfinx as df +import numpy as np +import ufl +from mpi4py import MPI + +from fenics_constitutive import SubSpaceMap, _build_view, build_subspace_map + + +def test_subspace_map_vector_equals_tensor(): + mesh = df.mesh.create_unit_cube(MPI.COMM_WORLD, 5, 7, 11) + + map_c = mesh.topology.index_map(mesh.topology.dim) + num_cells = map_c.size_local + map_c.num_ghosts + cells = np.arange(0, num_cells, dtype=np.int32) + + Q = ufl.FiniteElement("Quadrature", mesh.ufl_cell(), 2, quad_scheme="default") + QV = ufl.VectorElement("Quadrature", mesh.ufl_cell(), 2, quad_scheme="default", dim=3) + QT = ufl.TensorElement("Quadrature", mesh.ufl_cell(), 2, quad_scheme="default", shape=(3, 3)) + + Q_space = df.fem.FunctionSpace(mesh, Q) + QV_space = df.fem.FunctionSpace(mesh, QV) + QT_space = df.fem.FunctionSpace(mesh, QT) + + Q_map, _ = build_subspace_map(cells, Q_space) + QV_map, _ = _build_view(cells, QV_space) + QT_map, _ = _build_view(cells, QT_space) + + assert np.all(Q_map.parent==QV_map.parent) + assert np.all(Q_map.child==QV_map.child) + assert np.all(Q_map.parent==QT_map.parent) + assert np.all(Q_map.child==QT_map.child) + + for _ in range(10): + cell_sample = np.random.permutation(cells)[:num_cells//2] + + Q_map, _ = build_subspace_map(cells, Q_space) + QV_map, _ = build_subspace_map(cells, QV_space) + QT_map, _ = build_subspace_map(cells, QT_space) + + assert np.all(Q_map.parent==QV_map.parent) + assert np.all(Q_map.child==QV_map.child) + assert np.all(Q_map.parent==QT_map.parent) + assert np.all(Q_map.child==QT_map.child) + +def test_map_evaluation(): + mesh = df.mesh.create_unit_cube(MPI.COMM_WORLD, 5, 7, 11) + + map_c = mesh.topology.index_map(mesh.topology.dim) + num_cells = map_c.size_local + map_c.num_ghosts + cells = np.arange(0, num_cells, dtype=np.int32) + + Q = ufl.FiniteElement("Quadrature", mesh.ufl_cell(), 2, quad_scheme="default") + QV = ufl.VectorElement("Quadrature", mesh.ufl_cell(), 2, quad_scheme="default", dim=3) + QT = ufl.TensorElement("Quadrature", mesh.ufl_cell(), 2, quad_scheme="default", shape=(3, 3)) + + Q_space = df.fem.FunctionSpace(mesh, Q) + QV_space = df.fem.FunctionSpace(mesh, QV) + QT_space = df.fem.FunctionSpace(mesh, QT) + + q = df.fem.Function(Q_space) + q_test = q.copy() + qv = df.fem.Function(QV_space) + qv_test = qv.copy() + qt = df.fem.Function(QT_space) + qt_test = qt.copy() + + q.x.array[:] = np.random.random(q.x.array.shape) + qv.x.array[:] = np.random.random(qv.x.array.shape) + qt.x.array[:] = np.random.random(qt.x.array.shape) + + + for _ in range(10): + cell_sample = np.random.permutation(cells)[:num_cells//2] + + Q_sub_map, submesh = build_subspace_map(cell_sample, Q_space, return_subspace=False) + Q_sub = df.fem.FunctionSpace(submesh, Q) + QV_sub = df.fem.FunctionSpace(submesh, QV) + QT_sub = df.fem.FunctionSpace(submesh, QT) + + q_sub = df.fem.Function(Q_sub) + qv_sub = df.fem.Function(QV_sub) + qt_sub = df.fem.Function(QT_sub) + + Q_sub_map.map_to_child(q, q_sub) + Q_sub_map.map_to_child(qv, qv_sub) + Q_sub_map.map_to_child(qt, qt_sub) + + Q_sub_map.map_to_parent(q_sub, q_test) + Q_sub_map.map_to_parent(qv_sub, qv_test) + Q_sub_map.map_to_parent(qt_sub, qt_test) + + q_view = q.x.array.reshape(num_cells, -1) + q_test_view = q_test.x.array.reshape(num_cells, -1) + assert np.all(q_view[cell_sample] == q_test_view[cell_sample]) + + qv_view = qv.x.array.reshape(num_cells, -1) + qv_test_view = qv_test.x.array.reshape(num_cells, -1) + assert np.all(qv_view[cell_sample] == qv_test_view[cell_sample]) + + qt_view = qt.x.array.reshape(num_cells, -1) + qt_test_view = qt_test.x.array.reshape(num_cells, -1) + assert np.all(qt_view[cell_sample] == qt_test_view[cell_sample]) \ No newline at end of file From c762a752230cfcef4027e52e25c61dc856b82b3d Mon Sep 17 00:00:00 2001 From: srosenbu Date: Wed, 20 Mar 2024 10:17:46 +0100 Subject: [PATCH 06/19] include the examples in the pytest path --- pyproject.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/pyproject.toml b/pyproject.toml index 67d0ed9..63245ef 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -74,6 +74,7 @@ filterwarnings = [ log_cli_level = "INFO" testpaths = [ "tests", + "examples", ] From 78075d9ac2123d06c1a323cae32a9edc6368bfcd Mon Sep 17 00:00:00 2001 From: srosenbu Date: Wed, 20 Mar 2024 10:18:17 +0100 Subject: [PATCH 07/19] fix a bug --- src/fenics_constitutive/maps.py | 76 +++++---------------------------- 1 file changed, 11 insertions(+), 65 deletions(-) diff --git a/src/fenics_constitutive/maps.py b/src/fenics_constitutive/maps.py index 748f098..a7761b5 100644 --- a/src/fenics_constitutive/maps.py +++ b/src/fenics_constitutive/maps.py @@ -1,9 +1,11 @@ +from __future__ import annotations + +from dataclasses import dataclass + import dolfinx as df import numpy as np -from functools import reduce -from dataclasses import dataclass -__all__ = ["SubSpaceMap", "_build_view", "build_subspace_map"] +__all__ = ["SubSpaceMap", "build_subspace_map"] @dataclass @@ -19,12 +21,13 @@ def map_to_parent(self, sub: df.fem.Function, parent: df.fem.Function) -> None: parent.function_space.mesh == self.parent_mesh ), "Parent mesh does not match" assert sub.ufl_shape == parent.ufl_shape, "Shapes do not match" - # print(sub.ufl_element().value_size(), "sub.ufl_element().value_size()") - size = sub.ufl_element().value_size() # reduce(lambda x, y: x * y, sub.ufl_shape) if not isinstance(sub.ufl_shape, int) else sub.ufl_shape + + size = sub.ufl_element().value_size() parent_array = parent.x.array.reshape(-1, size) sub_array = sub.x.array.reshape(-1, size) parent_array[self.parent] = sub_array[self.child] + parent.x.scatter_forward() def map_to_child(self, parent: df.fem.Function, sub: df.fem.Function) -> None: assert sub.function_space.mesh == self.sub_mesh, "Subspace mesh does not match" @@ -33,12 +36,12 @@ def map_to_child(self, parent: df.fem.Function, sub: df.fem.Function) -> None: ), "Parent mesh does not match" assert sub.ufl_shape == parent.ufl_shape, "Shapes do not match" - # size = reduce(lambda x, y: x * y, sub.ufl_shape) if not isinstance(sub.ufl_shape, int) else sub.ufl_shape - size = sub.ufl_element().value_size() # reduce(lambda x, y: x * y, sub.ufl_shape) if not isinstance(sub.ufl_shape, int) else sub.ufl_shape + size = sub.ufl_element().value_size() parent_array = parent.x.array.reshape(-1, size) sub_array = sub.x.array.reshape(-1, size) - parent_array[self.parent] = sub_array[self.child] + sub_array[self.child] = parent_array[self.parent] + sub.x.scatter_forward() def build_subspace_map( @@ -80,60 +83,3 @@ def build_subspace_map( else: del V_sub return map, submesh - - -def _build_view( - cells: np.ndarray, V: df.fem.FunctionSpace -) -> tuple[SubSpaceMap, df.fem.FunctionSpace]: - mesh = V.mesh - submesh, cell_map, _, _ = df.mesh.create_submesh(mesh, mesh.topology.dim, cells) - fe = V.ufl_element() - V_sub = df.fem.FunctionSpace(submesh, fe) - - submesh = V_sub.mesh - view_parent = [] - view_child = [] - - num_sub_cells = submesh.topology.index_map(submesh.topology.dim).size_local - for cell in range(num_sub_cells): - view_child.append(V_sub.dofmap.cell_dofs(cell)) - view_parent.append(V.dofmap.cell_dofs(cell_map[cell])) - if view_child: - return ( - SubSpaceMap( - parent=np.hstack(view_parent), - child=np.hstack(view_child), - sub_mesh=submesh, - parent_mesh=V.mesh, - ), - V_sub, - ) - else: - # it may be that a process does not own any of the cells in the submesh - return ( - SubSpaceMap( - parent=np.array([], dtype=np.int32), - child=np.array([], dtype=np.int32), - sub_mesh=submesh, - parent_mesh=V.mesh, - ), - V_sub, - ) - return ( - np.array([], dtype=np.int32), - np.array([], dtype=np.int32), - V_sub, - ) - # if view_child: - # return ( - # np.hstack(view_parent), - # np.hstack(view_child), - # V_sub, - # ) - # else: - # # it may be that a process does not own any of the cells in the submesh - # return ( - # np.array([], dtype=np.int32), - # np.array([], dtype=np.int32), - # V_sub, - # ) From 43fc21d0d9e3fee1e3b79630c91ec61e50c83758 Mon Sep 17 00:00:00 2001 From: srosenbu Date: Wed, 20 Mar 2024 10:58:19 +0100 Subject: [PATCH 08/19] =?UTF-8?q?inc=C3=B6lude=20option=20for=20homogenous?= =?UTF-8?q?=20domains=20in=20solver?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/fenics_constitutive/solver.py | 200 ++++++++++++++++++++---------- tests/test_conversions.py | 22 ++-- tests/test_maps.py | 65 ++++++---- 3 files changed, 187 insertions(+), 100 deletions(-) diff --git a/src/fenics_constitutive/solver.py b/src/fenics_constitutive/solver.py index 9ba3205..d1b3687 100644 --- a/src/fenics_constitutive/solver.py +++ b/src/fenics_constitutive/solver.py @@ -1,28 +1,40 @@ -from dataclasses import dataclass -from typing import Any +from __future__ import annotations + import basix import dolfinx as df import numpy as np import ufl from petsc4py import PETSc + from .interfaces import IncrSmallStrainModel +from .maps import SubSpaceMap, build_subspace_map from .stress_strain import ufl_mandel_strain -from .maps import SubSpaceMap, build_subspace_map def build_history( - law: IncrSmallStrainModel, submesh: df.mesh.Mesh, q_degree: int + law: IncrSmallStrainModel, mesh: df.mesh.Mesh, q_degree: int ) -> df.fem.Function | dict[str, df.fem.Function] | None: + """Build the history space and function(s) for the given law. + + Args: + law: The constitutive law. + mesh: Either the full mesh for a homogenous domain or the submesh. + q_degree: The quadrature degree. + + Returns: + The history function(s) for the given law. + + """ match law.history_dim: case int(): Qh = ufl.VectorElement( "Quadrature", - submesh.ufl_cell(), + mesh.ufl_cell(), q_degree, quad_scheme="default", dim=law.history_dim, ) - history_space = df.fem.FunctionSpace(submesh, Qh) + history_space = df.fem.FunctionSpace(mesh, Qh) history = df.fem.Function(history_space) case None: history = None @@ -32,7 +44,7 @@ def build_history( if isinstance(value, int): Qh = ufl.VectorElement( "Quadrature", - submesh.ufl_cell(), + mesh.ufl_cell(), q_degree, quad_scheme="default", dim=value, @@ -40,41 +52,39 @@ def build_history( elif isinstance(value, tuple): Qh = ufl.TensorElement( "Quadrature", - submesh.ufl_cell(), + mesh.ufl_cell(), q_degree, quad_scheme="default", shape=value, ) - history_space = df.fem.FunctionSpace(submesh, Qh) + history_space = df.fem.FunctionSpace(mesh, Qh) history[key] = df.fem.Function(history_space) return history - -class IncrMechanicsProblem(df.fem.petsc.NonlinearProblem): +class IncrSmallStrainProblem(df.fem.petsc.NonlinearProblem): def __init__( self, - # laws = [(law_1, cells_1), (law_2, cells_2), ...] laws: list[tuple[IncrSmallStrainModel, np.ndarray]], - # u0 is basically a history variable - u0: df.fem.Function, + u: df.fem.Function, bcs: list[df.fem.DirichletBCMetaClass], q_degree: int = 1, - form_compiler_options: dict = {}, - jit_options: dict = {}, + form_compiler_options: dict | None = None, + jit_options: dict | None = None, ): - mesh = u0.function_space.mesh + mesh = u.function_space.mesh map_c = mesh.topology.index_map(mesh.topology.dim) num_cells = map_c.size_local + map_c.num_ghosts - self.cells = np.arange(0, num_cells, dtype=np.int32) - - self.gdim = mesh.ufl_cell().geometric_dimension() + cells = np.arange(0, num_cells, dtype=np.int32) constraint = laws[0][0].constraint assert all( law[0].constraint == constraint for law in laws ), "All laws must have the same constraint" + gdim = mesh.ufl_cell().geometric_dimension() + assert constraint.geometric_dim() == gdim, "Geometric dimension mismatch between mesh and laws" + QVe = ufl.VectorElement( "Quadrature", mesh.ufl_cell(), @@ -92,14 +102,10 @@ def __init__( QV = df.fem.FunctionSpace(mesh, QVe) QT = df.fem.FunctionSpace(mesh, QTe) - self.laws = [] - # TODO: We can probably reduce to one map per submesh - self.submesh_maps = [] - #self.QV_views = [] - #self.Qh_views = [] - #self.QT_views = [] + self.laws: list[tuple[IncrSmallStrainModel, np.ndarray]] = [] + self.submesh_maps: list[SubSpaceMap] = [] - self._strain = [] + self._del_strain = [] self._stress = [] self._history_0 = [] self._history_1 = [] @@ -108,22 +114,39 @@ def __init__( self._time = 0.0 # time at the end of the increment with df.common.Timer("submeshes-and-data-structures"): - for law, cells in laws: + if len(laws) > 1: + for law, cells in laws: + self.laws.append((law, cells)) + + # ### submesh and subspace for strain, stress + subspace_map, submesh, QV_subspace = build_subspace_map( + cells, QV, return_subspace=True + ) + self.submesh_maps.append(subspace_map) + self._del_strain.append(df.fem.Function(QV_subspace)) + self._stress.append(df.fem.Function(QV_subspace)) + + #subspace for tanget + QT_subspace = df.fem.FunctionSpace(submesh, QTe) + self._tangent.append(df.fem.Function(QT_subspace)) + + #subspaces for history + history_0 = build_history(law, submesh, q_degree) + history_1 = ( + {key: fn.copy() for key, fn in history_0.items()} + if isinstance(history_0, dict) + else history_0 + ) + self._history_0.append(history_0) + self._history_1.append(history_1) + else: + law, cells = laws[0] self.laws.append((law, cells)) + self._del_strain.append(df.fem.Function(QV)) - # ### submesh and subspace for strain, stress - subspace_map, submesh, QV_subspace = build_subspace_map(cells, QV, return_subspace=True) - self.submesh_maps.append(subspace_map) - self._strain.append(df.fem.Function(QV_subspace)) - self._stress.append(df.fem.Function(QV_subspace)) - - # ### submesh and subspace for tanget - QT_subspace = df.fem.FunctionSpace(submesh, QTe) - self._tangent.append(df.fem.Function(QT_subspace)) - - # ### submesh and subspace for history - history_0, Qh_views = build_history(law, submesh, q_degree) + #Spaces for history + history_0 = build_history(law, mesh, q_degree) history_1 = ( {key: fn.copy() for key, fn in history_0.items()} if isinstance(history_0, dict) @@ -136,15 +159,24 @@ def __init__( self.stress_1 = df.fem.Function(QV) self.tangent = df.fem.Function(QT) - u_, du = ufl.TestFunction(u0.function_space), ufl.TrialFunction(u0.function_space) + u_, du = ufl.TestFunction(u.function_space), ufl.TrialFunction(u.function_space) self.metadata = {"quadrature_degree": q_degree, "quadrature_scheme": "default"} self.dxm = ufl.dx(metadata=self.metadata) - self.R_form = ufl.inner(ufl_mandel_strain(u_, constraint), self.stress) * self.dxm - self.dR_form = ufl.inner(ufl_mandel_strain(du, constraint), ufl.dot(self.tangent, ufl_mandel_strain(u_, constraint))) * self.dxm + self.R_form = ( + ufl.inner(ufl_mandel_strain(u_, constraint), self.stress_1) * self.dxm + ) + self.dR_form = ( + ufl.inner( + ufl_mandel_strain(du, constraint), + ufl.dot(self.tangent, ufl_mandel_strain(u_, constraint)), + ) + * self.dxm + ) - self._u0 = u0 + self._u = u + self._u0 = u.copy() self._bcs = bcs self._form_compiler_options = form_compiler_options self._jit_options = jit_options @@ -152,6 +184,10 @@ def __init__( basix_celltype = getattr(basix.CellType, mesh.topology.cell_type.name) self.q_points, _ = basix.make_quadrature(basix_celltype, q_degree) + self.del_strain_expr = df.fem.Expression( + ufl_mandel_strain(self._u - self._u0, constraint), self.q_points + ) + # ### Note on JIT compilation of UFL forms # if super().__init__(R, u, bcs, dR) is called within ElasticityProblem.__init__ # the user cannot add Neumann BCs. @@ -172,8 +208,8 @@ def a(self) -> df.fem.FormMetaClass: self._u, self._bcs, self.dR_form, - form_compiler_options=self._form_compiler_options, - jit_options=self._jit_options, + form_compiler_options=self._form_compiler_options if self._form_compiler_options is not None else {}, + jit_options=self._jit_options if self._jit_options is not None else {}, ) return self._a @@ -187,34 +223,66 @@ def form(self, x: PETSc.Vec): The vector containing the latest solution """ super().form(x) - #TODO: look if the creation of the expression slows things down - strain_expr = df.fem.Expression(ufl_mandel_strain(x - self._u0, self.constraint), self.q_points) - for k, (law, cells) in enumerate(self.laws): + assert ( + x is self._u.vector + ), "The solution vector must be the same as the one passed to the MechanicsProblem" + if len(self.laws) > 1: + for k, (law, cells) in enumerate(self.laws): + with df.common.Timer("strain_evaluation"): + #TODO: test this!! + self.del_strain_expr.eval(cells, self._del_strain[k].x.array.reshape(cells.size, -1)) + + with df.common.Timer("stress_evaluation"): + self.submesh_maps[k].map_to_child(self.stress_0, self._stress[k]) + self._history_1[k].x.array[:] = self._history_0[k].x.array + law.evaluate( + self._time, + self._del_strain[k].x.array, + self._stress[k].x.array, + self._tangent[k].x.array, + self._history_1[k].x.array, # history, + ) + + with df.common.Timer("stress-local-to-global"): + self.submesh_maps[k].map_to_parent(self._stress[k], self.stress_1) + self.submesh_maps[k].map_to_parent(self._tangent[k], self.tangent) + else: + law, cells = self.laws[0] with df.common.Timer("strain_evaluation"): - strain_expr.eval( - cells, self._strain[k].x.array.reshape(cells.size, -1) - ) + self.del_strain_expr.eval(cells, self._del_strain[0].x.array) with df.common.Timer("stress_evaluation"): - self.submesh_maps[k].map_to_child(self.stress_0, self._stress[k]) - self._history_1[k].x.array[:] = self._history_0[k].x.array + self.stress_1.x.array[:] = self.stress_0.x.array + self._history_1[0].x.array[:] = self._history_0[0].x.array law.evaluate( self._time, - self._strain[k].x.array, - self._stress[k].x.array, - self._tangent[k].x.array, - self._history_1[k].x.array, # history, + self._del_strain[0].x.array, + self.stress_1.x.array, + self.tangent.x.array, + self._history_1[0].x.array, # history, ) - - with df.common.Timer("stress-local-to-global"): - self.submesh_maps[k].map_to_parent(self._stress[k], self.stress_1) - self.submesh_maps[k].map_to_parent(self._tangent[k], self.tangent) - - self.stress.x.scatter_forward() + self.stress_1.x.scatter_forward() self.tangent.x.scatter_forward() def update(self): - pass - + # update the current displacement, stress and history + # works for both homogeneous and inhomogeneous domains + self._u0.x.array[:] = self._u.x.array + self._u0.x.scatter_forward() + + self.stress_0.x.array[:] = self.stress_1.x.array + self.stress_0.x.scatter_forward() + + for k, (law, _) in enumerate(self.laws): + match law.history_dim: + case int(): + self._history_0[k].x.array[:] = self._history_1[k].x.array + self._history_0[k].x.scatter_forward() + case None: + pass + case dict(): + for key in law.history_dim: + self._history_0[k][key].x.array[:] = self._history_1[k][key].x.array + self._history_0[k][key].x.scatter_forward() diff --git a/tests/test_conversions.py b/tests/test_conversions.py index a6aa8de..d9bb6af 100644 --- a/tests/test_conversions.py +++ b/tests/test_conversions.py @@ -1,10 +1,13 @@ -from mpi4py import MPI -import ufl -from fenics_constitutive import Constraint, strain_from_grad_u, ufl_mandel_strain +from __future__ import annotations + +import basix +import dolfinx as df import numpy as np import pytest -import dolfinx as df -import basix +import ufl +from mpi4py import MPI + +from fenics_constitutive import Constraint, strain_from_grad_u, ufl_mandel_strain def test_strain_from_grad_u(): @@ -54,13 +57,16 @@ def test_ufl_strain_equals_array_conversion(constraint: Constraint): match constraint.geometric_dim(): case 1: mesh = df.mesh.create_unit_interval(MPI.COMM_WORLD, 2) - lam = lambda x: x[0] * 0.1 + def lam(x): + return x[0] * 0.1 case 2: mesh = df.mesh.create_unit_square(MPI.COMM_WORLD, 2, 2) - lam = lambda x: (x[0] * 0.1, x[1] * 0.2) + def lam(x): + return x[0] * 0.1, x[1] * 0.2 case 3: mesh = df.mesh.create_unit_cube(MPI.COMM_WORLD, 2, 2, 2) - lam = lambda x: (x[0] * 0.1, x[1] * 0.2, x[2] * 0.3) + def lam(x): + return x[0] * 0.1, x[1] * 0.2, x[2] * 0.3 P1 = df.fem.VectorFunctionSpace(mesh, ("Lagrange", 1)) u = df.fem.Function(P1) grad_u_ufl = ufl.grad(u) diff --git a/tests/test_maps.py b/tests/test_maps.py index 836bc76..e4a2c46 100644 --- a/tests/test_maps.py +++ b/tests/test_maps.py @@ -6,7 +6,7 @@ from fenics_constitutive import SubSpaceMap, _build_view, build_subspace_map -def test_subspace_map_vector_equals_tensor(): +def test_subspace_vector_map_vector_equals_tensor_map(): mesh = df.mesh.create_unit_cube(MPI.COMM_WORLD, 5, 7, 11) map_c = mesh.topology.index_map(mesh.topology.dim) @@ -14,8 +14,12 @@ def test_subspace_map_vector_equals_tensor(): cells = np.arange(0, num_cells, dtype=np.int32) Q = ufl.FiniteElement("Quadrature", mesh.ufl_cell(), 2, quad_scheme="default") - QV = ufl.VectorElement("Quadrature", mesh.ufl_cell(), 2, quad_scheme="default", dim=3) - QT = ufl.TensorElement("Quadrature", mesh.ufl_cell(), 2, quad_scheme="default", shape=(3, 3)) + QV = ufl.VectorElement( + "Quadrature", mesh.ufl_cell(), 2, quad_scheme="default", dim=3 + ) + QT = ufl.TensorElement( + "Quadrature", mesh.ufl_cell(), 2, quad_scheme="default", shape=(3, 3) + ) Q_space = df.fem.FunctionSpace(mesh, Q) QV_space = df.fem.FunctionSpace(mesh, QV) @@ -25,22 +29,23 @@ def test_subspace_map_vector_equals_tensor(): QV_map, _ = _build_view(cells, QV_space) QT_map, _ = _build_view(cells, QT_space) - assert np.all(Q_map.parent==QV_map.parent) - assert np.all(Q_map.child==QV_map.child) - assert np.all(Q_map.parent==QT_map.parent) - assert np.all(Q_map.child==QT_map.child) - + assert np.all(Q_map.parent == QV_map.parent) + assert np.all(Q_map.child == QV_map.child) + assert np.all(Q_map.parent == QT_map.parent) + assert np.all(Q_map.child == QT_map.child) + for _ in range(10): - cell_sample = np.random.permutation(cells)[:num_cells//2] + cell_sample = np.random.permutation(cells)[: num_cells // 2] Q_map, _ = build_subspace_map(cells, Q_space) QV_map, _ = build_subspace_map(cells, QV_space) QT_map, _ = build_subspace_map(cells, QT_space) - assert np.all(Q_map.parent==QV_map.parent) - assert np.all(Q_map.child==QV_map.child) - assert np.all(Q_map.parent==QT_map.parent) - assert np.all(Q_map.child==QT_map.child) + assert np.all(Q_map.parent == QV_map.parent) + assert np.all(Q_map.child == QV_map.child) + assert np.all(Q_map.parent == QT_map.parent) + assert np.all(Q_map.child == QT_map.child) + def test_map_evaluation(): mesh = df.mesh.create_unit_cube(MPI.COMM_WORLD, 5, 7, 11) @@ -50,8 +55,12 @@ def test_map_evaluation(): cells = np.arange(0, num_cells, dtype=np.int32) Q = ufl.FiniteElement("Quadrature", mesh.ufl_cell(), 2, quad_scheme="default") - QV = ufl.VectorElement("Quadrature", mesh.ufl_cell(), 2, quad_scheme="default", dim=3) - QT = ufl.TensorElement("Quadrature", mesh.ufl_cell(), 2, quad_scheme="default", shape=(3, 3)) + QV = ufl.VectorElement( + "Quadrature", mesh.ufl_cell(), 2, quad_scheme="default", dim=3 + ) + QT = ufl.TensorElement( + "Quadrature", mesh.ufl_cell(), 2, quad_scheme="default", shape=(3, 3) + ) Q_space = df.fem.FunctionSpace(mesh, Q) QV_space = df.fem.FunctionSpace(mesh, QV) @@ -64,15 +73,19 @@ def test_map_evaluation(): qt = df.fem.Function(QT_space) qt_test = qt.copy() - q.x.array[:] = np.random.random(q.x.array.shape) - qv.x.array[:] = np.random.random(qv.x.array.shape) - qt.x.array[:] = np.random.random(qt.x.array.shape) - + scalar_array = np.random.random(q.x.array.shape) + q.x.array[:] = scalar_array + vector_array = np.random.random(qv.x.array.shape) + qv.x.array[:] = vector_array + tensor_array = np.random.random(qt.x.array.shape) + qt.x.array[:] = tensor_array for _ in range(10): - cell_sample = np.random.permutation(cells)[:num_cells//2] + cell_sample = np.random.permutation(cells)[: num_cells // 2] - Q_sub_map, submesh = build_subspace_map(cell_sample, Q_space, return_subspace=False) + Q_sub_map, submesh = build_subspace_map( + cell_sample, Q_space, return_subspace=False + ) Q_sub = df.fem.FunctionSpace(submesh, Q) QV_sub = df.fem.FunctionSpace(submesh, QV) QT_sub = df.fem.FunctionSpace(submesh, QT) @@ -89,14 +102,14 @@ def test_map_evaluation(): Q_sub_map.map_to_parent(qv_sub, qv_test) Q_sub_map.map_to_parent(qt_sub, qt_test) - q_view = q.x.array.reshape(num_cells, -1) + q_view = scalar_array.reshape(num_cells, -1) q_test_view = q_test.x.array.reshape(num_cells, -1) assert np.all(q_view[cell_sample] == q_test_view[cell_sample]) - - qv_view = qv.x.array.reshape(num_cells, -1) + + qv_view = vector_array.reshape(num_cells, -1) qv_test_view = qv_test.x.array.reshape(num_cells, -1) assert np.all(qv_view[cell_sample] == qv_test_view[cell_sample]) - qt_view = qt.x.array.reshape(num_cells, -1) + qt_view = tensor_array.reshape(num_cells, -1) qt_test_view = qt_test.x.array.reshape(num_cells, -1) - assert np.all(qt_view[cell_sample] == qt_test_view[cell_sample]) \ No newline at end of file + assert np.all(qt_view[cell_sample] == qt_test_view[cell_sample]) From fb6fcc54e7b483df1d63a6e6ed3f0e5424a85127 Mon Sep 17 00:00:00 2001 From: srosenbu Date: Wed, 20 Mar 2024 13:18:57 +0100 Subject: [PATCH 09/19] First test for solver and linear elasticity --- .../linear_elasticity_model.py | 14 +++--- examples/linear_elasticity/test_elasticity.py | 49 +++++++++++++++++++ src/fenics_constitutive/solver.py | 41 +++++++++++----- src/fenics_constitutive/stress_strain.py | 5 +- tests/test_maps.py | 6 +-- 5 files changed, 91 insertions(+), 24 deletions(-) create mode 100644 examples/linear_elasticity/test_elasticity.py diff --git a/examples/linear_elasticity/linear_elasticity_model.py b/examples/linear_elasticity/linear_elasticity_model.py index 4bbb00e..694ce40 100644 --- a/examples/linear_elasticity/linear_elasticity_model.py +++ b/examples/linear_elasticity/linear_elasticity_model.py @@ -69,15 +69,15 @@ def evaluate( history: np.ndarray | dict[str, np.ndarray] | None, ) -> None: assert ( - grad_del_u.size / (self.geometric_dim**2) - == mandel_stress.size / self.stress_strain_dim - == tangent.size / (self.stress_strain_dim**2) + grad_del_u.size // (self.geometric_dim**2) + == mandel_stress.size // self.stress_strain_dim + == tangent.size // (self.stress_strain_dim**2) ) - n_gauss = grad_del_u.shape / (self.geometric_dim**2) - + n_gauss = grad_del_u.size // (self.geometric_dim**2) + mandel_view = mandel_stress.reshape(-1, self.stress_strain_dim) strain_increment = strain_from_grad_u(grad_del_u, self.constraint) - mandel_stress += strain_increment.reshape(-1, self.stress_strain_dim) @ self.D - tangent[:] = np.tile(self.D, n_gauss) + mandel_view += strain_increment.reshape(-1, self.stress_strain_dim) @ self.D + tangent[:] = np.tile(self.D.flatten(), n_gauss) @property def constraint(self) -> Constraint: diff --git a/examples/linear_elasticity/test_elasticity.py b/examples/linear_elasticity/test_elasticity.py new file mode 100644 index 0000000..56c620e --- /dev/null +++ b/examples/linear_elasticity/test_elasticity.py @@ -0,0 +1,49 @@ +from mpi4py import MPI +import dolfinx as df +from linear_elasticity_model import LinearElasticityModel +from fenics_constitutive import Constraint, IncrSmallStrainProblem +import numpy as np +from dolfinx.nls.petsc import NewtonSolver + +youngs_modulus = 42. +poissons_ratio = 0.3 + +def test_uniaxial_stress(): + mesh = df.mesh.create_unit_interval(MPI.COMM_WORLD, 2) + V = df.fem.FunctionSpace(mesh, ("CG", 1)) + u = df.fem.Function(V) + law = LinearElasticityModel( + parameters={"E": youngs_modulus, "nu": poissons_ratio}, + constraint=Constraint.UNIAXIAL_STRESS, + ) + + def left_boundary(x): + return np.isclose(x[0], 0.0) + + def right_boundary(x): + return np.isclose(x[0], 1.0) + + displacement = df.fem.Constant(mesh, 0.01) + dofs_left = df.fem.locate_dofs_geometrical(V, left_boundary) + dofs_right = df.fem.locate_dofs_geometrical(V, right_boundary) + bc_left = df.fem.dirichletbc(df.fem.Constant(mesh,0.0), dofs_left, V) + bc_right = df.fem.dirichletbc(displacement, dofs_right, V) + + + problem = IncrSmallStrainProblem( + law, + u, + [bc_left, bc_right], + ) + + solver = NewtonSolver(MPI.COMM_WORLD, problem) + n, converged = solver.solve(u) + assert abs(problem.stress_1.x.array[0] - youngs_modulus * 0.01) < 1e-10 / (youngs_modulus * 0.01) + + problem.update() + assert abs(problem.stress_0.x.array[0] - youngs_modulus * 0.01) < 1e-10 / (youngs_modulus * 0.01) + assert np.max(problem._u0.x.array) == displacement.value + + displacement.value = 0.02 + n, converged = solver.solve(u) + assert abs(problem.stress_1.x.array[0] - youngs_modulus * 0.02) < 1e-10 / (youngs_modulus * 0.02) diff --git a/src/fenics_constitutive/solver.py b/src/fenics_constitutive/solver.py index d1b3687..ea3f44c 100644 --- a/src/fenics_constitutive/solver.py +++ b/src/fenics_constitutive/solver.py @@ -65,7 +65,7 @@ def build_history( class IncrSmallStrainProblem(df.fem.petsc.NonlinearProblem): def __init__( self, - laws: list[tuple[IncrSmallStrainModel, np.ndarray]], + laws: list[tuple[IncrSmallStrainModel, np.ndarray]] | IncrSmallStrainModel, u: df.fem.Function, bcs: list[df.fem.DirichletBCMetaClass], q_degree: int = 1, @@ -76,6 +76,8 @@ def __init__( map_c = mesh.topology.index_map(mesh.topology.dim) num_cells = map_c.size_local + map_c.num_ghosts cells = np.arange(0, num_cells, dtype=np.int32) + if isinstance(laws, IncrSmallStrainModel): + laws = [(laws, cells)] constraint = laws[0][0].constraint assert all( @@ -99,13 +101,20 @@ def __init__( quad_scheme="default", shape=(constraint.stress_strain_dim(), constraint.stress_strain_dim()), ) + Q_grad_u_e = ufl.TensorElement( + "Quadrature", + mesh.ufl_cell(), + q_degree, + quad_scheme="default", + shape=(gdim, gdim), + ) QV = df.fem.FunctionSpace(mesh, QVe) QT = df.fem.FunctionSpace(mesh, QTe) self.laws: list[tuple[IncrSmallStrainModel, np.ndarray]] = [] self.submesh_maps: list[SubSpaceMap] = [] - self._del_strain = [] + self._del_grad_u = [] self._stress = [] self._history_0 = [] self._history_1 = [] @@ -123,8 +132,11 @@ def __init__( cells, QV, return_subspace=True ) self.submesh_maps.append(subspace_map) - self._del_strain.append(df.fem.Function(QV_subspace)) self._stress.append(df.fem.Function(QV_subspace)) + + #subspace for grad u + Q_grad_u_subspace = df.fem.FunctionSpace(submesh, Q_grad_u_e) + self._del_grad_u.append(df.fem.Function(Q_grad_u_subspace)) #subspace for tanget QT_subspace = df.fem.FunctionSpace(submesh, QTe) @@ -143,7 +155,9 @@ def __init__( law, cells = laws[0] self.laws.append((law, cells)) - self._del_strain.append(df.fem.Function(QV)) + #subspace for grad u + Q_grad_u_space = df.fem.FunctionSpace(mesh, Q_grad_u_e) + self._del_grad_u.append(df.fem.Function(Q_grad_u_space)) #Spaces for history history_0 = build_history(law, mesh, q_degree) @@ -184,8 +198,8 @@ def __init__( basix_celltype = getattr(basix.CellType, mesh.topology.cell_type.name) self.q_points, _ = basix.make_quadrature(basix_celltype, q_degree) - self.del_strain_expr = df.fem.Expression( - ufl_mandel_strain(self._u - self._u0, constraint), self.q_points + self.del_grad_u_expr = df.fem.Expression( + ufl.nabla_grad(self._u - self._u0), self.q_points ) # ### Note on JIT compilation of UFL forms @@ -225,20 +239,20 @@ def form(self, x: PETSc.Vec): super().form(x) assert ( - x is self._u.vector + x.array.data == self._u.vector.array.data ), "The solution vector must be the same as the one passed to the MechanicsProblem" if len(self.laws) > 1: for k, (law, cells) in enumerate(self.laws): with df.common.Timer("strain_evaluation"): #TODO: test this!! - self.del_strain_expr.eval(cells, self._del_strain[k].x.array.reshape(cells.size, -1)) + self.del_grad_u_expr.eval(cells, self._del_grad_u[k].x.array.reshape(cells.size, -1)) with df.common.Timer("stress_evaluation"): self.submesh_maps[k].map_to_child(self.stress_0, self._stress[k]) self._history_1[k].x.array[:] = self._history_0[k].x.array law.evaluate( self._time, - self._del_strain[k].x.array, + self._del_grad_u[k].x.array, self._stress[k].x.array, self._tangent[k].x.array, self._history_1[k].x.array, # history, @@ -250,17 +264,18 @@ def form(self, x: PETSc.Vec): else: law, cells = self.laws[0] with df.common.Timer("strain_evaluation"): - self.del_strain_expr.eval(cells, self._del_strain[0].x.array) + self.del_grad_u_expr.eval(cells, self._del_grad_u[0].x.array.reshape(cells.size, -1)) with df.common.Timer("stress_evaluation"): self.stress_1.x.array[:] = self.stress_0.x.array - self._history_1[0].x.array[:] = self._history_0[0].x.array + if law.history_dim is not None: + self._history_1[0].x.array[:] = self._history_0[0].x.array law.evaluate( self._time, - self._del_strain[0].x.array, + self._del_grad_u[0].x.array, self.stress_1.x.array, self.tangent.x.array, - self._history_1[0].x.array, # history, + self._history_1[0].x.array if law.history_dim is not None else None, # history, ) self.stress_1.x.scatter_forward() diff --git a/src/fenics_constitutive/stress_strain.py b/src/fenics_constitutive/stress_strain.py index c168607..02f2e15 100644 --- a/src/fenics_constitutive/stress_strain.py +++ b/src/fenics_constitutive/stress_strain.py @@ -20,7 +20,10 @@ def ufl_mandel_strain( Vector-valued UFL expression of the mandel strain. """ strain_dim = constraint.stress_strain_dim() - assert u.ufl_shape == (constraint.geometric_dim(),) + + shape = len(u.ufl_shape) + geometric_dim = u.ufl_shape[0] if shape > 0 else 1 + assert geometric_dim == constraint.geometric_dim() match constraint: case Constraint.UNIAXIAL_STRAIN: return ufl.nabla_grad(u) diff --git a/tests/test_maps.py b/tests/test_maps.py index e4a2c46..44fb31e 100644 --- a/tests/test_maps.py +++ b/tests/test_maps.py @@ -3,7 +3,7 @@ import ufl from mpi4py import MPI -from fenics_constitutive import SubSpaceMap, _build_view, build_subspace_map +from fenics_constitutive import SubSpaceMap, build_subspace_map def test_subspace_vector_map_vector_equals_tensor_map(): @@ -26,8 +26,8 @@ def test_subspace_vector_map_vector_equals_tensor_map(): QT_space = df.fem.FunctionSpace(mesh, QT) Q_map, _ = build_subspace_map(cells, Q_space) - QV_map, _ = _build_view(cells, QV_space) - QT_map, _ = _build_view(cells, QT_space) + QV_map, _ = build_subspace_map(cells, QV_space) + QT_map, _ = build_subspace_map(cells, QT_space) assert np.all(Q_map.parent == QV_map.parent) assert np.all(Q_map.child == QV_map.child) From 4f1e9bf5646f837584eb99b0a67bcee50b422b1f Mon Sep 17 00:00:00 2001 From: srosenbu Date: Wed, 20 Mar 2024 13:53:38 +0100 Subject: [PATCH 10/19] Test for two different constitutive laws --- examples/linear_elasticity/test_elasticity.py | 82 ++++++++++++++++--- src/fenics_constitutive/solver.py | 5 +- 2 files changed, 74 insertions(+), 13 deletions(-) diff --git a/examples/linear_elasticity/test_elasticity.py b/examples/linear_elasticity/test_elasticity.py index 56c620e..f8b5d7b 100644 --- a/examples/linear_elasticity/test_elasticity.py +++ b/examples/linear_elasticity/test_elasticity.py @@ -1,13 +1,17 @@ -from mpi4py import MPI +from __future__ import annotations + import dolfinx as df -from linear_elasticity_model import LinearElasticityModel -from fenics_constitutive import Constraint, IncrSmallStrainProblem import numpy as np from dolfinx.nls.petsc import NewtonSolver +from linear_elasticity_model import LinearElasticityModel +from mpi4py import MPI -youngs_modulus = 42. +from fenics_constitutive import Constraint, IncrSmallStrainProblem + +youngs_modulus = 42.0 poissons_ratio = 0.3 + def test_uniaxial_stress(): mesh = df.mesh.create_unit_interval(MPI.COMM_WORLD, 2) V = df.fem.FunctionSpace(mesh, ("CG", 1)) @@ -19,17 +23,16 @@ def test_uniaxial_stress(): def left_boundary(x): return np.isclose(x[0], 0.0) - + def right_boundary(x): return np.isclose(x[0], 1.0) - + displacement = df.fem.Constant(mesh, 0.01) dofs_left = df.fem.locate_dofs_geometrical(V, left_boundary) dofs_right = df.fem.locate_dofs_geometrical(V, right_boundary) - bc_left = df.fem.dirichletbc(df.fem.Constant(mesh,0.0), dofs_left, V) + bc_left = df.fem.dirichletbc(df.fem.Constant(mesh, 0.0), dofs_left, V) bc_right = df.fem.dirichletbc(displacement, dofs_right, V) - problem = IncrSmallStrainProblem( law, u, @@ -38,12 +41,69 @@ def right_boundary(x): solver = NewtonSolver(MPI.COMM_WORLD, problem) n, converged = solver.solve(u) - assert abs(problem.stress_1.x.array[0] - youngs_modulus * 0.01) < 1e-10 / (youngs_modulus * 0.01) + assert abs(problem.stress_1.x.array[0] - youngs_modulus * 0.01) < 1e-10 / ( + youngs_modulus * 0.01 + ) problem.update() - assert abs(problem.stress_0.x.array[0] - youngs_modulus * 0.01) < 1e-10 / (youngs_modulus * 0.01) + assert abs(problem.stress_0.x.array[0] - youngs_modulus * 0.01) < 1e-10 / ( + youngs_modulus * 0.01 + ) assert np.max(problem._u0.x.array) == displacement.value displacement.value = 0.02 n, converged = solver.solve(u) - assert abs(problem.stress_1.x.array[0] - youngs_modulus * 0.02) < 1e-10 / (youngs_modulus * 0.02) + assert abs(problem.stress_1.x.array[0] - youngs_modulus * 0.02) < 1e-10 / ( + youngs_modulus * 0.02 + ) + + +def test_uniaxial_stress_two_laws(): + mesh = df.mesh.create_unit_interval(MPI.COMM_WORLD, 2) + V = df.fem.FunctionSpace(mesh, ("CG", 1)) + u = df.fem.Function(V) + laws = [ + ( + LinearElasticityModel( + parameters={"E": youngs_modulus, "nu": poissons_ratio}, + constraint=Constraint.UNIAXIAL_STRESS, + ), + np.array([0], dtype=np.int32), + ), + ( + LinearElasticityModel( + parameters={"E": 2.0 * youngs_modulus, "nu": poissons_ratio}, + constraint=Constraint.UNIAXIAL_STRESS, + ), + np.array([1], dtype=np.int32), + ), + ] + + def left_boundary(x): + return np.isclose(x[0], 0.0) + + def right_boundary(x): + return np.isclose(x[0], 1.0) + + displacement = df.fem.Constant(mesh, 0.01) + dofs_left = df.fem.locate_dofs_geometrical(V, left_boundary) + dofs_right = df.fem.locate_dofs_geometrical(V, right_boundary) + bc_left = df.fem.dirichletbc(df.fem.Constant(mesh, 0.0), dofs_left, V) + bc_right = df.fem.dirichletbc(displacement, dofs_right, V) + + problem = IncrSmallStrainProblem( + laws, + u, + [bc_left, bc_right], + ) + + solver = NewtonSolver(MPI.COMM_WORLD, problem) + n, converged = solver.solve(u) + problem.update() + assert abs(problem.stress_0.x.array[0] - problem.stress_0.x.array[1]) < 1e-10 / abs( + problem.stress_0.x.array[0] + ) + + assert abs( + problem._del_grad_u[0].x.array[0] - 2.0 * problem._del_grad_u[1].x.array[0] + ) < 1e-10 / abs(problem._del_grad_u[0].x.array[0]) diff --git a/src/fenics_constitutive/solver.py b/src/fenics_constitutive/solver.py index ea3f44c..27de696 100644 --- a/src/fenics_constitutive/solver.py +++ b/src/fenics_constitutive/solver.py @@ -249,13 +249,14 @@ def form(self, x: PETSc.Vec): with df.common.Timer("stress_evaluation"): self.submesh_maps[k].map_to_child(self.stress_0, self._stress[k]) - self._history_1[k].x.array[:] = self._history_0[k].x.array + if law.history_dim is not None: + self._history_1[0].x.array[:] = self._history_0[0].x.array law.evaluate( self._time, self._del_grad_u[k].x.array, self._stress[k].x.array, self._tangent[k].x.array, - self._history_1[k].x.array, # history, + self._history_1[k].x.array if law.history_dim is not None else None, ) with df.common.Timer("stress-local-to-global"): From 7eb7407c37a329d44643438d069b68f8e1c9c230 Mon Sep 17 00:00:00 2001 From: srosenbu Date: Wed, 20 Mar 2024 13:58:59 +0100 Subject: [PATCH 11/19] More test cases --- examples/linear_elasticity/test_elasticity.py | 20 +++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/examples/linear_elasticity/test_elasticity.py b/examples/linear_elasticity/test_elasticity.py index f8b5d7b..bb4726c 100644 --- a/examples/linear_elasticity/test_elasticity.py +++ b/examples/linear_elasticity/test_elasticity.py @@ -2,6 +2,7 @@ import dolfinx as df import numpy as np +import pytest from dolfinx.nls.petsc import NewtonSolver from linear_elasticity_model import LinearElasticityModel from mpi4py import MPI @@ -57,8 +58,16 @@ def right_boundary(x): youngs_modulus * 0.02 ) - -def test_uniaxial_stress_two_laws(): +@pytest.mark.parametrize( + ["factor"], + [ + [0.5], + [2.0], + [3.0], + [4.0], + ], +) +def test_uniaxial_stress_two_laws(factor: float): mesh = df.mesh.create_unit_interval(MPI.COMM_WORLD, 2) V = df.fem.FunctionSpace(mesh, ("CG", 1)) u = df.fem.Function(V) @@ -72,7 +81,7 @@ def test_uniaxial_stress_two_laws(): ), ( LinearElasticityModel( - parameters={"E": 2.0 * youngs_modulus, "nu": poissons_ratio}, + parameters={"E": factor * youngs_modulus, "nu": poissons_ratio}, constraint=Constraint.UNIAXIAL_STRESS, ), np.array([1], dtype=np.int32), @@ -100,10 +109,13 @@ def right_boundary(x): solver = NewtonSolver(MPI.COMM_WORLD, problem) n, converged = solver.solve(u) problem.update() + + # Is the stress homogenous? assert abs(problem.stress_0.x.array[0] - problem.stress_0.x.array[1]) < 1e-10 / abs( problem.stress_0.x.array[0] ) + # Does the stiffer element have a proportionally lower strain? assert abs( - problem._del_grad_u[0].x.array[0] - 2.0 * problem._del_grad_u[1].x.array[0] + problem._del_grad_u[0].x.array[0] - factor * problem._del_grad_u[1].x.array[0] ) < 1e-10 / abs(problem._del_grad_u[0].x.array[0]) From db647ba6bf7151c72c1992da5ec65ec54f57f646 Mon Sep 17 00:00:00 2001 From: srosenbu Date: Wed, 20 Mar 2024 14:07:56 +0100 Subject: [PATCH 12/19] Uniaxial strain test --- examples/linear_elasticity/test_elasticity.py | 43 +++++++++++++++++++ 1 file changed, 43 insertions(+) diff --git a/examples/linear_elasticity/test_elasticity.py b/examples/linear_elasticity/test_elasticity.py index bb4726c..4dc15ef 100644 --- a/examples/linear_elasticity/test_elasticity.py +++ b/examples/linear_elasticity/test_elasticity.py @@ -58,6 +58,7 @@ def right_boundary(x): youngs_modulus * 0.02 ) + @pytest.mark.parametrize( ["factor"], [ @@ -119,3 +120,45 @@ def right_boundary(x): assert abs( problem._del_grad_u[0].x.array[0] - factor * problem._del_grad_u[1].x.array[0] ) < 1e-10 / abs(problem._del_grad_u[0].x.array[0]) + + +def test_uniaxial_strain(): + mesh = df.mesh.create_unit_interval(MPI.COMM_WORLD, 2) + V = df.fem.FunctionSpace(mesh, ("CG", 1)) + u = df.fem.Function(V) + law = LinearElasticityModel( + parameters={"E": youngs_modulus, "nu": poissons_ratio}, + constraint=Constraint.UNIAXIAL_STRAIN, + ) + + def left_boundary(x): + return np.isclose(x[0], 0.0) + + def right_boundary(x): + return np.isclose(x[0], 1.0) + + displacement = df.fem.Constant(mesh, 0.01) + dofs_left = df.fem.locate_dofs_geometrical(V, left_boundary) + dofs_right = df.fem.locate_dofs_geometrical(V, right_boundary) + bc_left = df.fem.dirichletbc(df.fem.Constant(mesh, 0.0), dofs_left, V) + bc_right = df.fem.dirichletbc(displacement, dofs_right, V) + + problem = IncrSmallStrainProblem( + law, + u, + [bc_left, bc_right], + ) + + solver = NewtonSolver(MPI.COMM_WORLD, problem) + n, converged = solver.solve(u) + problem.update() + + analytical_stress = ( + youngs_modulus + * (1.0 - poissons_ratio) + / ((1.0 + poissons_ratio) * (1.0 - 2.0 * poissons_ratio)) + ) * displacement.value + + assert abs(problem.stress_0.x.array[0] - analytical_stress) < 1e-10 / ( + analytical_stress + ) From ba6d2792877fdeaa1fa97ff25bec929ed008323a Mon Sep 17 00:00:00 2001 From: srosenbu Date: Wed, 20 Mar 2024 15:15:18 +0100 Subject: [PATCH 13/19] 3d test case --- examples/linear_elasticity/test_elasticity.py | 59 +++++++++++++++++++ src/fenics_constitutive/solver.py | 50 ++++++++++------ 2 files changed, 91 insertions(+), 18 deletions(-) diff --git a/examples/linear_elasticity/test_elasticity.py b/examples/linear_elasticity/test_elasticity.py index 4dc15ef..35c7e18 100644 --- a/examples/linear_elasticity/test_elasticity.py +++ b/examples/linear_elasticity/test_elasticity.py @@ -3,6 +3,7 @@ import dolfinx as df import numpy as np import pytest +import ufl from dolfinx.nls.petsc import NewtonSolver from linear_elasticity_model import LinearElasticityModel from mpi4py import MPI @@ -162,3 +163,61 @@ def right_boundary(x): assert abs(problem.stress_0.x.array[0] - analytical_stress) < 1e-10 / ( analytical_stress ) + + +def test_3d(): + # test the 3d case against a pure fenics implementation + mesh = df.mesh.create_unit_cube(MPI.COMM_WORLD, 2, 2, 2) + V = df.fem.VectorFunctionSpace(mesh, ("CG", 1)) + u = df.fem.Function(V) + law = LinearElasticityModel( + parameters={"E": youngs_modulus, "nu": poissons_ratio}, + constraint=Constraint.FULL, + ) + + def left_boundary(x): + return np.isclose(x[0], 0.0) + + def right_boundary(x): + return np.isclose(x[0], 1.0) + + # displacement = df.fem.Constant(mesh_2d, np.ndarray([0.01,0.0])) + dofs_left = df.fem.locate_dofs_geometrical(V, left_boundary) + dofs_right = df.fem.locate_dofs_geometrical(V, right_boundary) + bc_left = df.fem.dirichletbc(np.array([0.0,0.0, 0.0]), dofs_left, V) + bc_right = df.fem.dirichletbc(np.array([0.01, 0.0, 0.0]), dofs_right, V) + + problem = IncrSmallStrainProblem( + law, + u, + [bc_left, bc_right], + ) + + solver = NewtonSolver(MPI.COMM_WORLD, problem) + n, converged = solver.solve(u) + problem.update() + + v_, u_ = ufl.TestFunction(V), ufl.TrialFunction(V) + + def eps(v): + return ufl.sym(ufl.grad(v)) + + def sigma(v): + eps = ufl.sym(ufl.grad(v)) + lam, mu = ( + youngs_modulus + * poissons_ratio + / ((1.0 + poissons_ratio) * (1.0 - 2.0 * poissons_ratio)), + youngs_modulus / (2.0 * (1.0 + poissons_ratio)), + ) + return 2.0 * mu * eps + lam * ufl.tr(eps) * ufl.Identity(len(v)) + zero = df.fem.Constant(mesh, np.array([0.0, 0.0, 0.0])) + a = ufl.inner(ufl.grad(v_), sigma(u_)) * ufl.dx + L = ufl.dot(zero ,v_) * ufl.dx + + u_fenics = u.copy() + problem_fenics = df.fem.petsc.LinearProblem(a, L, [bc_left, bc_right], u=u_fenics) + problem_fenics.solve() + + # Check that the solution is the same + assert np.linalg.norm(u_fenics.x.array - u.x.array) < 1e-8/np.linalg.norm(u_fenics.x.array) \ No newline at end of file diff --git a/src/fenics_constitutive/solver.py b/src/fenics_constitutive/solver.py index 27de696..7716a91 100644 --- a/src/fenics_constitutive/solver.py +++ b/src/fenics_constitutive/solver.py @@ -15,12 +15,12 @@ def build_history( law: IncrSmallStrainModel, mesh: df.mesh.Mesh, q_degree: int ) -> df.fem.Function | dict[str, df.fem.Function] | None: """Build the history space and function(s) for the given law. - + Args: law: The constitutive law. mesh: Either the full mesh for a homogenous domain or the submesh. q_degree: The quadrature degree. - + Returns: The history function(s) for the given law. @@ -85,7 +85,9 @@ def __init__( ), "All laws must have the same constraint" gdim = mesh.ufl_cell().geometric_dimension() - assert constraint.geometric_dim() == gdim, "Geometric dimension mismatch between mesh and laws" + assert ( + constraint.geometric_dim() == gdim + ), "Geometric dimension mismatch between mesh and laws" QVe = ufl.VectorElement( "Quadrature", @@ -133,16 +135,16 @@ def __init__( ) self.submesh_maps.append(subspace_map) self._stress.append(df.fem.Function(QV_subspace)) - - #subspace for grad u + + # subspace for grad u Q_grad_u_subspace = df.fem.FunctionSpace(submesh, Q_grad_u_e) self._del_grad_u.append(df.fem.Function(Q_grad_u_subspace)) - #subspace for tanget + # subspace for tanget QT_subspace = df.fem.FunctionSpace(submesh, QTe) self._tangent.append(df.fem.Function(QT_subspace)) - #subspaces for history + # subspaces for history history_0 = build_history(law, submesh, q_degree) history_1 = ( {key: fn.copy() for key, fn in history_0.items()} @@ -154,12 +156,12 @@ def __init__( else: law, cells = laws[0] self.laws.append((law, cells)) - - #subspace for grad u + + # subspace for grad u Q_grad_u_space = df.fem.FunctionSpace(mesh, Q_grad_u_e) self._del_grad_u.append(df.fem.Function(Q_grad_u_space)) - #Spaces for history + # Spaces for history history_0 = build_history(law, mesh, q_degree) history_1 = ( {key: fn.copy() for key, fn in history_0.items()} @@ -222,7 +224,9 @@ def a(self) -> df.fem.FormMetaClass: self._u, self._bcs, self.dR_form, - form_compiler_options=self._form_compiler_options if self._form_compiler_options is not None else {}, + form_compiler_options=self._form_compiler_options + if self._form_compiler_options is not None + else {}, jit_options=self._jit_options if self._jit_options is not None else {}, ) @@ -244,8 +248,10 @@ def form(self, x: PETSc.Vec): if len(self.laws) > 1: for k, (law, cells) in enumerate(self.laws): with df.common.Timer("strain_evaluation"): - #TODO: test this!! - self.del_grad_u_expr.eval(cells, self._del_grad_u[k].x.array.reshape(cells.size, -1)) + # TODO: test this!! + self.del_grad_u_expr.eval( + cells, self._del_grad_u[k].x.array.reshape(cells.size, -1) + ) with df.common.Timer("stress_evaluation"): self.submesh_maps[k].map_to_child(self.stress_0, self._stress[k]) @@ -256,7 +262,9 @@ def form(self, x: PETSc.Vec): self._del_grad_u[k].x.array, self._stress[k].x.array, self._tangent[k].x.array, - self._history_1[k].x.array if law.history_dim is not None else None, + self._history_1[k].x.array + if law.history_dim is not None + else None, ) with df.common.Timer("stress-local-to-global"): @@ -265,7 +273,9 @@ def form(self, x: PETSc.Vec): else: law, cells = self.laws[0] with df.common.Timer("strain_evaluation"): - self.del_grad_u_expr.eval(cells, self._del_grad_u[0].x.array.reshape(cells.size, -1)) + self.del_grad_u_expr.eval( + cells, self._del_grad_u[0].x.array.reshape(cells.size, -1) + ) with df.common.Timer("stress_evaluation"): self.stress_1.x.array[:] = self.stress_0.x.array @@ -276,7 +286,9 @@ def form(self, x: PETSc.Vec): self._del_grad_u[0].x.array, self.stress_1.x.array, self.tangent.x.array, - self._history_1[0].x.array if law.history_dim is not None else None, # history, + self._history_1[0].x.array + if law.history_dim is not None + else None, # history, ) self.stress_1.x.scatter_forward() @@ -290,7 +302,7 @@ def update(self): self.stress_0.x.array[:] = self.stress_1.x.array self.stress_0.x.scatter_forward() - + for k, (law, _) in enumerate(self.laws): match law.history_dim: case int(): @@ -300,5 +312,7 @@ def update(self): pass case dict(): for key in law.history_dim: - self._history_0[k][key].x.array[:] = self._history_1[k][key].x.array + self._history_0[k][key].x.array[:] = self._history_1[k][ + key + ].x.array self._history_0[k][key].x.scatter_forward() From a1a4c5e8b2b1a5d3d6fbc1f124086fe834b3e639 Mon Sep 17 00:00:00 2001 From: srosenbu Date: Thu, 21 Mar 2024 09:21:12 +0100 Subject: [PATCH 14/19] Fix an error in the linear elasticity model --- examples/linear_elasticity/linear_elasticity_model.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/linear_elasticity/linear_elasticity_model.py b/examples/linear_elasticity/linear_elasticity_model.py index 694ce40..605ab7b 100644 --- a/examples/linear_elasticity/linear_elasticity_model.py +++ b/examples/linear_elasticity/linear_elasticity_model.py @@ -27,8 +27,8 @@ def __init__(self, parameters: dict[str, float], constraint: Constraint): # see https://en.wikipedia.org/wiki/Hooke%27s_law self.D = np.array( [ - [2.0 * mu + lam, lam, 0.0, 0.0], - [lam, 2.0 * mu + lam, 0.0, 0.0], + [2.0 * mu + lam, lam, lam, 0.0], + [lam, 2.0 * mu + lam, lam, 0.0], [lam, lam, 2.0 * mu + lam, 0.0], [0.0, 0.0, 0.0, 2.0 * mu], ] From dd9e188cc14ec2622cad1b0d9f91d0aee486f5c2 Mon Sep 17 00:00:00 2001 From: srosenbu Date: Thu, 21 Mar 2024 09:22:02 +0100 Subject: [PATCH 15/19] Sanity checks for all constraints --- examples/linear_elasticity/test_elasticity.py | 90 ++++++++++++++++++- 1 file changed, 87 insertions(+), 3 deletions(-) diff --git a/examples/linear_elasticity/test_elasticity.py b/examples/linear_elasticity/test_elasticity.py index 35c7e18..df99b2e 100644 --- a/examples/linear_elasticity/test_elasticity.py +++ b/examples/linear_elasticity/test_elasticity.py @@ -164,6 +164,87 @@ def right_boundary(x): analytical_stress ) +def test_plane_strain(): + # sanity check if out of plane stress is NOT zero + mesh = df.mesh.create_unit_square(MPI.COMM_WORLD, 2, 2) + V = df.fem.VectorFunctionSpace(mesh, ("CG", 1)) + u = df.fem.Function(V) + law = LinearElasticityModel( + parameters={"E": youngs_modulus, "nu": poissons_ratio}, + constraint=Constraint.PLANE_STRAIN, + ) + + def left_boundary(x): + return np.isclose(x[0], 0.0) + + def right_boundary(x): + return np.isclose(x[0], 1.0) + + dofs_left = df.fem.locate_dofs_geometrical(V, left_boundary) + dofs_right = df.fem.locate_dofs_geometrical(V, right_boundary) + bc_left = df.fem.dirichletbc(np.array([0.0, 0.0]), dofs_left, V) + bc_right = df.fem.dirichletbc(np.array([0.01, 0.0]), dofs_right, V) + + problem = IncrSmallStrainProblem( + law, + u, + [bc_left, bc_right], + ) + + solver = NewtonSolver(MPI.COMM_WORLD, problem) + n, converged = solver.solve(u) + problem.update() + print(problem.stress_0.x.array.reshape(-1, law.constraint.stress_strain_dim())) + assert ( + np.linalg.norm( + problem.stress_0.x.array.reshape(-1, law.constraint.stress_strain_dim())[ + :, 2 + ] + ) + > 1e-2 + ) + + +def test_plane_stress(): + # just a sanity check if out of plane stress is really zero + mesh = df.mesh.create_unit_square(MPI.COMM_WORLD, 2, 2) + V = df.fem.VectorFunctionSpace(mesh, ("CG", 1)) + u = df.fem.Function(V) + law = LinearElasticityModel( + parameters={"E": youngs_modulus, "nu": poissons_ratio}, + constraint=Constraint.PLANE_STRESS, + ) + + def left_boundary(x): + return np.isclose(x[0], 0.0) + + def right_boundary(x): + return np.isclose(x[0], 1.0) + + # displacement = df.fem.Constant(mesh_2d, np.ndarray([0.01,0.0])) + dofs_left = df.fem.locate_dofs_geometrical(V, left_boundary) + dofs_right = df.fem.locate_dofs_geometrical(V, right_boundary) + bc_left = df.fem.dirichletbc(np.array([0.0, 0.0]), dofs_left, V) + bc_right = df.fem.dirichletbc(np.array([0.01, 0.0]), dofs_right, V) + + problem = IncrSmallStrainProblem( + law, + u, + [bc_left, bc_right], + ) + + solver = NewtonSolver(MPI.COMM_WORLD, problem) + n, converged = solver.solve(u) + problem.update() + assert ( + np.linalg.norm( + problem.stress_0.x.array.reshape(-1, law.constraint.stress_strain_dim())[ + :, 2 + ] + ) + < 1e-10 + ) + def test_3d(): # test the 3d case against a pure fenics implementation @@ -184,7 +265,7 @@ def right_boundary(x): # displacement = df.fem.Constant(mesh_2d, np.ndarray([0.01,0.0])) dofs_left = df.fem.locate_dofs_geometrical(V, left_boundary) dofs_right = df.fem.locate_dofs_geometrical(V, right_boundary) - bc_left = df.fem.dirichletbc(np.array([0.0,0.0, 0.0]), dofs_left, V) + bc_left = df.fem.dirichletbc(np.array([0.0, 0.0, 0.0]), dofs_left, V) bc_right = df.fem.dirichletbc(np.array([0.01, 0.0, 0.0]), dofs_right, V) problem = IncrSmallStrainProblem( @@ -211,13 +292,16 @@ def sigma(v): youngs_modulus / (2.0 * (1.0 + poissons_ratio)), ) return 2.0 * mu * eps + lam * ufl.tr(eps) * ufl.Identity(len(v)) + zero = df.fem.Constant(mesh, np.array([0.0, 0.0, 0.0])) a = ufl.inner(ufl.grad(v_), sigma(u_)) * ufl.dx - L = ufl.dot(zero ,v_) * ufl.dx + L = ufl.dot(zero, v_) * ufl.dx u_fenics = u.copy() problem_fenics = df.fem.petsc.LinearProblem(a, L, [bc_left, bc_right], u=u_fenics) problem_fenics.solve() # Check that the solution is the same - assert np.linalg.norm(u_fenics.x.array - u.x.array) < 1e-8/np.linalg.norm(u_fenics.x.array) \ No newline at end of file + assert np.linalg.norm(u_fenics.x.array - u.x.array) < 1e-8 / np.linalg.norm( + u_fenics.x.array + ) From 005a34566fa0b1d6aa5de2113846a72890e0c817 Mon Sep 17 00:00:00 2001 From: srosenbu Date: Thu, 21 Mar 2024 09:39:19 +0100 Subject: [PATCH 16/19] Make ruff happy --- .../linear_elasticity_model.py | 10 +++-- examples/linear_elasticity/test_elasticity.py | 14 +++--- pyproject.toml | 2 +- src/fenics_constitutive/__init__.py | 4 +- src/fenics_constitutive/conversion.py | 44 ++++++++++--------- src/fenics_constitutive/interfaces.py | 6 +-- src/fenics_constitutive/stress_strain.py | 7 +-- tests/test_conversions.py | 26 +++++------ tests/test_maps.py | 10 +++-- 9 files changed, 66 insertions(+), 57 deletions(-) diff --git a/examples/linear_elasticity/linear_elasticity_model.py b/examples/linear_elasticity/linear_elasticity_model.py index 605ab7b..961b981 100644 --- a/examples/linear_elasticity/linear_elasticity_model.py +++ b/examples/linear_elasticity/linear_elasticity_model.py @@ -1,6 +1,9 @@ -from fenics_constitutive import Constraint, IncrSmallStrainModel, strain_from_grad_u +from __future__ import annotations + import numpy as np +from fenics_constitutive import Constraint, IncrSmallStrainModel, strain_from_grad_u + class LinearElasticityModel(IncrSmallStrainModel): def __init__(self, parameters: dict[str, float], constraint: Constraint): @@ -34,7 +37,7 @@ def __init__(self, parameters: dict[str, float], constraint: Constraint): ] ) case Constraint.PLANE_STRESS: - # We do not make any assumptions about strain in the z-direction + # We do not make any assumptions about strain in the z-direction # This matrix just multiplies the z component by 0.0 which results # in a plane stress state # see https://en.wikipedia.org/wiki/Hooke%27s_law @@ -58,7 +61,8 @@ def __init__(self, parameters: dict[str, float], constraint: Constraint): # see https://csmbrannon.net/2012/08/02/distinction-between-uniaxial-stress-and-uniaxial-strain/ self.D = np.array([[E]]) case _: - raise NotImplementedError("Constraint not implemented") + msg = "Constraint not implemented" + raise NotImplementedError(msg) def evaluate( self, diff --git a/examples/linear_elasticity/test_elasticity.py b/examples/linear_elasticity/test_elasticity.py index df99b2e..ffc4e8e 100644 --- a/examples/linear_elasticity/test_elasticity.py +++ b/examples/linear_elasticity/test_elasticity.py @@ -61,12 +61,12 @@ def right_boundary(x): @pytest.mark.parametrize( - ["factor"], + ("factor"), [ - [0.5], - [2.0], - [3.0], - [4.0], + (0.5), + (2.0), + (3.0), + (4.0), ], ) def test_uniaxial_stress_two_laws(factor: float): @@ -164,6 +164,7 @@ def right_boundary(x): analytical_stress ) + def test_plane_strain(): # sanity check if out of plane stress is NOT zero mesh = df.mesh.create_unit_square(MPI.COMM_WORLD, 2, 2) @@ -179,7 +180,7 @@ def left_boundary(x): def right_boundary(x): return np.isclose(x[0], 1.0) - + dofs_left = df.fem.locate_dofs_geometrical(V, left_boundary) dofs_right = df.fem.locate_dofs_geometrical(V, right_boundary) bc_left = df.fem.dirichletbc(np.array([0.0, 0.0]), dofs_left, V) @@ -194,7 +195,6 @@ def right_boundary(x): solver = NewtonSolver(MPI.COMM_WORLD, problem) n, converged = solver.solve(u) problem.update() - print(problem.stress_0.x.array.reshape(-1, law.constraint.stress_strain_dim())) assert ( np.linalg.norm( problem.stress_0.x.array.reshape(-1, law.constraint.stress_strain_dim())[ diff --git a/pyproject.toml b/pyproject.toml index 63245ef..9dde59e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -103,7 +103,7 @@ disallow_incomplete_defs = true [tool.ruff] -src = ["src"] +src = ["src", "tests", "examples"] [tool.ruff.lint] extend-select = [ diff --git a/src/fenics_constitutive/__init__.py b/src/fenics_constitutive/__init__.py index f001bb6..2dc7c6b 100644 --- a/src/fenics_constitutive/__init__.py +++ b/src/fenics_constitutive/__init__.py @@ -9,8 +9,8 @@ from ._version import version as __version__ from .interfaces import * -from .stress_strain import * -from .solver import * from .maps import * +from .solver import * +from .stress_strain import * __all__ = ["__version__"] diff --git a/src/fenics_constitutive/conversion.py b/src/fenics_constitutive/conversion.py index 6161917..e84e06e 100644 --- a/src/fenics_constitutive/conversion.py +++ b/src/fenics_constitutive/conversion.py @@ -1,8 +1,12 @@ +from __future__ import annotations + import numpy as np + from .interfaces import Constraint __all__ = ["strain_from_grad_u"] + def strain_from_grad_u(grad_u: np.ndarray, constraint: Constraint) -> np.ndarray: """Compute the Mandel-strain from the gradient of displacement. @@ -13,16 +17,16 @@ def strain_from_grad_u(grad_u: np.ndarray, constraint: Constraint) -> np.ndarray Strain tensor. """ strain_dim = constraint.stress_strain_dim() - n_gauss = int(grad_u.size / (constraint.geometric_dim()**2)) - strain = np.zeros(strain_dim*n_gauss) - grad_u_view = grad_u.reshape(-1, constraint.geometric_dim()**2) + n_gauss = int(grad_u.size / (constraint.geometric_dim() ** 2)) + strain = np.zeros(strain_dim * n_gauss) + grad_u_view = grad_u.reshape(-1, constraint.geometric_dim() ** 2) strain_view = strain.reshape(-1, strain_dim) match constraint: case Constraint.UNIAXIAL_STRAIN: strain_view[:, 0] = grad_u_view[:, 0] case Constraint.UNIAXIAL_STRESS: - strain_view[:,0] = grad_u_view[:, 0] + strain_view[:, 0] = grad_u_view[:, 0] case Constraint.PLANE_STRAIN: """ Full tensor: @@ -34,10 +38,10 @@ def strain_from_grad_u(grad_u: np.ndarray, constraint: Constraint) -> np.ndarray f = 1 / sqrt(2) 0 3 "0" f*(1+2) """ - strain_view[:,0] = grad_u_view[:,0] - strain_view[:,1] = grad_u_view[:,3] - strain_view[:,2] = 0.0 - strain_view[:,3] = 1/2**0.5 * (grad_u_view[:,1] + grad_u_view[:,2]) + strain_view[:, 0] = grad_u_view[:, 0] + strain_view[:, 1] = grad_u_view[:, 3] + strain_view[:, 2] = 0.0 + strain_view[:, 3] = 1 / 2**0.5 * (grad_u_view[:, 1] + grad_u_view[:, 2]) case Constraint.PLANE_STRESS: """ Full tensor: @@ -49,10 +53,10 @@ def strain_from_grad_u(grad_u: np.ndarray, constraint: Constraint) -> np.ndarray f = 1 / sqrt(2) 0 3 "0" f*(1+2) """ - strain_view[:,0] = grad_u_view[:,0] - strain_view[:,1] = grad_u_view[:,3] - strain_view[:,2] = 0.0 - strain_view[:,3] = 1/2**0.5 * (grad_u_view[:,1] + grad_u_view[:,2]) + strain_view[:, 0] = grad_u_view[:, 0] + strain_view[:, 1] = grad_u_view[:, 3] + strain_view[:, 2] = 0.0 + strain_view[:, 3] = 1 / 2**0.5 * (grad_u_view[:, 1] + grad_u_view[:, 2]) case Constraint.FULL: """ Full tensor: @@ -65,13 +69,13 @@ def strain_from_grad_u(grad_u: np.ndarray, constraint: Constraint) -> np.ndarray f = 1 / sqrt(2) 0 4 8 f*(1+3) f*(5+7) f*(2+6) """ - strain_view[:,0] = grad_u_view[:,0] - strain_view[:,1] = grad_u_view[:,4] - strain_view[:,2] = grad_u_view[:,8] - strain_view[:,3] = 1/2**0.5 * (grad_u_view[:,1] + grad_u_view[:,3]) - strain_view[:,4] = 1/2**0.5 * (grad_u_view[:,5] + grad_u_view[:,7]) - strain_view[:,5] = 1/2**0.5 * (grad_u_view[:,2] + grad_u_view[:,6]) + strain_view[:, 0] = grad_u_view[:, 0] + strain_view[:, 1] = grad_u_view[:, 4] + strain_view[:, 2] = grad_u_view[:, 8] + strain_view[:, 3] = 1 / 2**0.5 * (grad_u_view[:, 1] + grad_u_view[:, 3]) + strain_view[:, 4] = 1 / 2**0.5 * (grad_u_view[:, 5] + grad_u_view[:, 7]) + strain_view[:, 5] = 1 / 2**0.5 * (grad_u_view[:, 2] + grad_u_view[:, 6]) case _: - raise NotImplementedError("Only full constraint implemented") + msg = "Only full constraint implemented" + raise NotImplementedError(msg) return strain - \ No newline at end of file diff --git a/src/fenics_constitutive/interfaces.py b/src/fenics_constitutive/interfaces.py index 453b1f6..bdceccb 100644 --- a/src/fenics_constitutive/interfaces.py +++ b/src/fenics_constitutive/interfaces.py @@ -1,3 +1,5 @@ +from __future__ import annotations + from abc import ABC, abstractmethod, abstractproperty from enum import Enum @@ -90,7 +92,6 @@ def evaluate( tangent : The tangent. history : The history variable(s). """ - pass @abstractmethod def update(self) -> None: @@ -100,7 +101,6 @@ def update(self) -> None: For example: The model could contain the current time which is not stored in the history, but needs to be updated after each evaluation. """ - pass @abstractproperty def constraint(self) -> Constraint: @@ -110,7 +110,6 @@ def constraint(self) -> Constraint: Returns: The constraint. """ - pass @property def stress_strain_dim(self) -> int: @@ -144,4 +143,3 @@ def history_dim(self) -> int | dict[str, int | tuple[int, int]] | None: Returns: The dimension of the history variable(s). """ - pass diff --git a/src/fenics_constitutive/stress_strain.py b/src/fenics_constitutive/stress_strain.py index 02f2e15..38301d0 100644 --- a/src/fenics_constitutive/stress_strain.py +++ b/src/fenics_constitutive/stress_strain.py @@ -1,3 +1,5 @@ +from __future__ import annotations + import numpy as np import ufl @@ -19,8 +21,6 @@ def ufl_mandel_strain( Returns: Vector-valued UFL expression of the mandel strain. """ - strain_dim = constraint.stress_strain_dim() - shape = len(u.ufl_shape) geometric_dim = u.ufl_shape[0] if shape > 0 else 1 assert geometric_dim == constraint.geometric_dim() @@ -135,5 +135,6 @@ def strain_from_grad_u(grad_u: np.ndarray, constraint: Constraint) -> np.ndarray strain_view[:, 4] = 1 / 2**0.5 * (grad_u_view[:, 5] + grad_u_view[:, 7]) strain_view[:, 5] = 1 / 2**0.5 * (grad_u_view[:, 2] + grad_u_view[:, 6]) case _: - raise NotImplementedError("Constraint not supported.") + msg = "Constraint not supported." + raise NotImplementedError(msg) return strain diff --git a/tests/test_conversions.py b/tests/test_conversions.py index d9bb6af..94b54a0 100644 --- a/tests/test_conversions.py +++ b/tests/test_conversions.py @@ -44,29 +44,33 @@ def test_strain_from_grad_u(): @pytest.mark.parametrize( - ["constraint"], + ("constraint"), [ - [Constraint.UNIAXIAL_STRAIN], - [Constraint.UNIAXIAL_STRESS], - [Constraint.PLANE_STRAIN], - [Constraint.PLANE_STRESS], - [Constraint.FULL], + (Constraint.UNIAXIAL_STRAIN), + (Constraint.UNIAXIAL_STRESS), + (Constraint.PLANE_STRAIN), + (Constraint.PLANE_STRESS), + (Constraint.FULL), ], ) def test_ufl_strain_equals_array_conversion(constraint: Constraint): match constraint.geometric_dim(): case 1: mesh = df.mesh.create_unit_interval(MPI.COMM_WORLD, 2) + def lam(x): return x[0] * 0.1 case 2: mesh = df.mesh.create_unit_square(MPI.COMM_WORLD, 2, 2) + def lam(x): return x[0] * 0.1, x[1] * 0.2 case 3: mesh = df.mesh.create_unit_cube(MPI.COMM_WORLD, 2, 2, 2) + def lam(x): return x[0] * 0.1, x[1] * 0.2, x[2] * 0.3 + P1 = df.fem.VectorFunctionSpace(mesh, ("Lagrange", 1)) u = df.fem.Function(P1) grad_u_ufl = ufl.grad(u) @@ -79,13 +83,9 @@ def lam(x): 1, ) - expr_grad_u = df.fem.Expression( - grad_u_ufl, points - ) - expr_mandel_strain = df.fem.Expression( - mandel_strain_ufl, points - ) - + expr_grad_u = df.fem.Expression(grad_u_ufl, points) + expr_mandel_strain = df.fem.Expression(mandel_strain_ufl, points) + n_cells = mesh.topology.index_map(mesh.topology.dim).size_local grad_u_array = expr_grad_u.eval(np.arange(n_cells, dtype=np.int32)).flatten() strain_array = expr_mandel_strain.eval(np.arange(n_cells, dtype=np.int32)).flatten() diff --git a/tests/test_maps.py b/tests/test_maps.py index 44fb31e..0ab0c4e 100644 --- a/tests/test_maps.py +++ b/tests/test_maps.py @@ -1,9 +1,11 @@ +from __future__ import annotations + import dolfinx as df import numpy as np import ufl from mpi4py import MPI -from fenics_constitutive import SubSpaceMap, build_subspace_map +from fenics_constitutive import build_subspace_map def test_subspace_vector_map_vector_equals_tensor_map(): @@ -37,9 +39,9 @@ def test_subspace_vector_map_vector_equals_tensor_map(): for _ in range(10): cell_sample = np.random.permutation(cells)[: num_cells // 2] - Q_map, _ = build_subspace_map(cells, Q_space) - QV_map, _ = build_subspace_map(cells, QV_space) - QT_map, _ = build_subspace_map(cells, QT_space) + Q_map, _ = build_subspace_map(cell_sample, Q_space) + QV_map, _ = build_subspace_map(cell_sample, QV_space) + QT_map, _ = build_subspace_map(cell_sample, QT_space) assert np.all(Q_map.parent == QV_map.parent) assert np.all(Q_map.child == QV_map.child) From 429e1faf3f107a7be8397b46ad09989ffcb3b7ac Mon Sep 17 00:00:00 2001 From: srosenbu Date: Fri, 22 Mar 2024 09:21:26 +0100 Subject: [PATCH 17/19] docstrings --- environment.yml | 1 + examples/linear_elasticity/test_elasticity.py | 6 ++ src/fenics_constitutive/interfaces.py | 3 +- src/fenics_constitutive/maps.py | 37 +++++++++++++ src/fenics_constitutive/solver.py | 55 ++++++++++++------- 5 files changed, 80 insertions(+), 22 deletions(-) diff --git a/environment.yml b/environment.yml index 2564a99..3654617 100644 --- a/environment.yml +++ b/environment.yml @@ -8,5 +8,6 @@ dependencies: - numpy - pip - pytest + - pytest-mpi - pip: - -e . \ No newline at end of file diff --git a/examples/linear_elasticity/test_elasticity.py b/examples/linear_elasticity/test_elasticity.py index ffc4e8e..c570e06 100644 --- a/examples/linear_elasticity/test_elasticity.py +++ b/examples/linear_elasticity/test_elasticity.py @@ -43,18 +43,24 @@ def right_boundary(x): solver = NewtonSolver(MPI.COMM_WORLD, problem) n, converged = solver.solve(u) + + # Compare the result with the analytical solution assert abs(problem.stress_1.x.array[0] - youngs_modulus * 0.01) < 1e-10 / ( youngs_modulus * 0.01 ) problem.update() + # Check that the stress is updated correctly assert abs(problem.stress_0.x.array[0] - youngs_modulus * 0.01) < 1e-10 / ( youngs_modulus * 0.01 ) + # Check that the displacement is updated correctly assert np.max(problem._u0.x.array) == displacement.value displacement.value = 0.02 n, converged = solver.solve(u) + + # Compare the result of the updated problem with new BC with the analytical solution assert abs(problem.stress_1.x.array[0] - youngs_modulus * 0.02) < 1e-10 / ( youngs_modulus * 0.02 ) diff --git a/src/fenics_constitutive/interfaces.py b/src/fenics_constitutive/interfaces.py index bdceccb..7a2c550 100644 --- a/src/fenics_constitutive/interfaces.py +++ b/src/fenics_constitutive/interfaces.py @@ -83,7 +83,8 @@ def evaluate( tangent: np.ndarray, history: np.ndarray | dict[str, np.ndarray] | None, ) -> None: - """Evaluate the constitutive model and overwrite the stress, tangent and history. + """ + Evaluate the constitutive model and overwrite the stress, tangent and history. Args: del_t : The time increment. diff --git a/src/fenics_constitutive/maps.py b/src/fenics_constitutive/maps.py index a7761b5..b9420a8 100644 --- a/src/fenics_constitutive/maps.py +++ b/src/fenics_constitutive/maps.py @@ -10,12 +10,30 @@ @dataclass class SubSpaceMap: + """ + Map between a subspace and a parent space. + + Args: + parent: The dof indices of the parent space. + child: The dof indices of the child space. + sub_mesh: The mesh of the subspace. + parent_mesh: The mesh of the parent space. + + """ + parent: np.ndarray child: np.ndarray sub_mesh: df.mesh.Mesh parent_mesh: df.mesh.Mesh def map_to_parent(self, sub: df.fem.Function, parent: df.fem.Function) -> None: + """ + Map the values from the subspace to the parent space. + + Args: + sub: The function in the subspace. + parent: The function in the parent space. + """ assert sub.function_space.mesh == self.sub_mesh, "Subspace mesh does not match" assert ( parent.function_space.mesh == self.parent_mesh @@ -30,6 +48,13 @@ def map_to_parent(self, sub: df.fem.Function, parent: df.fem.Function) -> None: parent.x.scatter_forward() def map_to_child(self, parent: df.fem.Function, sub: df.fem.Function) -> None: + """ + Map the values from the parent space to the subspace. + + Args: + parent: The function in the parent space. + sub: The function in the subspace. + """ assert sub.function_space.mesh == self.sub_mesh, "Subspace mesh does not match" assert ( parent.function_space.mesh == self.parent_mesh @@ -50,6 +75,18 @@ def build_subspace_map( tuple[SubSpaceMap, df.mesh.Mesh] | tuple[SubSpaceMap, df.mesh.Mesh, df.fem.FunctionSpace] ): + """ + Build a map between a subspace and a parent space. This currently needs + to build a functionspace which can optionally be returned. + + Args: + cells: The cells of the subspace. + V: The parent function space. + return_subspace: Return the subspace function space. + + Returns: + The subspace map, the submesh and optionally the subspace function space. + """ mesh = V.mesh submesh, cell_map, _, _ = df.mesh.create_submesh(mesh, mesh.topology.dim, cells) fe = V.ufl_element() diff --git a/src/fenics_constitutive/solver.py b/src/fenics_constitutive/solver.py index 7716a91..4e7b263 100644 --- a/src/fenics_constitutive/solver.py +++ b/src/fenics_constitutive/solver.py @@ -63,6 +63,28 @@ def build_history( class IncrSmallStrainProblem(df.fem.petsc.NonlinearProblem): + """ + A nonlinear problem for incremental small strain models. To be used with + the dolfinx NewtonSolver. + + Args: + laws: A list of tuples where the first element is the constitutive law and the second + element is the cells for the submesh. If only one law is provided, it is assumed + that the domain is homogenous. + u: The displacement field. This is the unknown in the nonlinear problem. + bcs: The Dirichlet boundary conditions. + q_degree: The quadrature degree (Polynomial degree which the quadrature rule needs to integrate exactly). + form_compiler_options: The options for the form compiler. + jit_options: The options for the JIT compiler. + + Note: + If `super().__init__(R, u, bcs, dR)` is called within the __init__ method, + the user cannot add Neumann BCs. Therefore, the compilation (i.e. call to + `super().__init__()`) is done when `df.nls.petsc.NewtonSolver` is initialized. + The solver will call `self._A = fem.petsc.create_matrix(problem.a)` and hence + we override the property ``a`` of NonlinearProblem to ensure that the form is compiled. + """ + def __init__( self, laws: list[tuple[IncrSmallStrainModel, np.ndarray]] | IncrSmallStrainModel, @@ -204,15 +226,6 @@ def __init__( ufl.nabla_grad(self._u - self._u0), self.q_points ) - # ### Note on JIT compilation of UFL forms - # if super().__init__(R, u, bcs, dR) is called within ElasticityProblem.__init__ - # the user cannot add Neumann BCs. - # Therefore, the compilation (i.e. call to super().__init__()) is done when - # df.nls.petsc.NewtonSolver is initialized. - # df.nls.petsc.NewtonSolver will call - # self._A = fem.petsc.create_matrix(problem.a) and hence it would suffice - # to override the property ``a`` of NonlinearProblem. - @property def a(self) -> df.fem.FormMetaClass: """Compiled bilinear form (the Jacobian form)""" @@ -232,13 +245,14 @@ def a(self) -> df.fem.FormMetaClass: return self._a - def form(self, x: PETSc.Vec): + def form(self, x: PETSc.Vec) -> None: """This function is called before the residual or Jacobian is - computed. This is usually used to update ghost values. - Parameters - ---------- - x - The vector containing the latest solution + computed. This is usually used to update ghost values, but here + we use it to update the stress, tangent and history. + + Args: + x: The vector containing the latest solution + """ super().form(x) @@ -286,17 +300,16 @@ def form(self, x: PETSc.Vec): self._del_grad_u[0].x.array, self.stress_1.x.array, self.tangent.x.array, - self._history_1[0].x.array - if law.history_dim is not None - else None, # history, + self._history_1[0].x.array if law.history_dim is not None else None, ) self.stress_1.x.scatter_forward() self.tangent.x.scatter_forward() - def update(self): - # update the current displacement, stress and history - # works for both homogeneous and inhomogeneous domains + def update(self) -> None: + """ + Update the current displacement, stress and history. + """ self._u0.x.array[:] = self._u.x.array self._u0.x.scatter_forward() From 309d95358a5e8c96deb1144e5d8604220d7c4c4f Mon Sep 17 00:00:00 2001 From: srosenbu Date: Fri, 22 Mar 2024 09:31:02 +0100 Subject: [PATCH 18/19] try mpi pytest --- .github/workflows/pytest.yml | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml index 838c3fd..f99f046 100644 --- a/.github/workflows/pytest.yml +++ b/.github/workflows/pytest.yml @@ -54,12 +54,11 @@ jobs: shell: bash -l {0} run: | pytest - #- name: execute example scripts - # shell: bash -l {0} - # run: | - # for file in docs/examples/*.py - # do - # python $file - # done + + - name: pytest-mpi + id : pytest-mpi + shell: bash -l {0} + run: | + mpirun -np 2 python3 -m pytest From 36ff2599d84e0cdfc528e749ffacb6e9e7225edb Mon Sep 17 00:00:00 2001 From: Philipp Diercks Date: Fri, 22 Mar 2024 12:20:38 +0100 Subject: [PATCH 19/19] Revert "try mpi pytest" This reverts commit 309d95358a5e8c96deb1144e5d8604220d7c4c4f. --- .github/workflows/pytest.yml | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml index f99f046..838c3fd 100644 --- a/.github/workflows/pytest.yml +++ b/.github/workflows/pytest.yml @@ -54,11 +54,12 @@ jobs: shell: bash -l {0} run: | pytest - - - name: pytest-mpi - id : pytest-mpi - shell: bash -l {0} - run: | - mpirun -np 2 python3 -m pytest + #- name: execute example scripts + # shell: bash -l {0} + # run: | + # for file in docs/examples/*.py + # do + # python $file + # done