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/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/linear_elasticity_model.py b/examples/linear_elasticity/linear_elasticity_model.py new file mode 100644 index 0000000..961b981 --- /dev/null +++ b/examples/linear_elasticity/linear_elasticity_model.py @@ -0,0 +1,95 @@ +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): + self._constraint = constraint + 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 + 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, 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, 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], + ] + ) + case Constraint.PLANE_STRESS: + # 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 + / (1 - nu**2.0) + * np.array( + [ + [1.0, nu, 0.0, 0.0], + [nu, 1.0, 0.0, 0.0], + [0.0, 0.0, 0.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 _: + msg = "Constraint not implemented" + raise NotImplementedError(msg) + + 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.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_view += strain_increment.reshape(-1, self.stress_strain_dim) @ self.D + tangent[:] = np.tile(self.D.flatten(), n_gauss) + + @property + def constraint(self) -> Constraint: + return self._constraint + + @property + def history_dim(self) -> None: + return None + + def update(self) -> None: + pass diff --git a/examples/linear_elasticity/test_elasticity.py b/examples/linear_elasticity/test_elasticity.py new file mode 100644 index 0000000..c570e06 --- /dev/null +++ b/examples/linear_elasticity/test_elasticity.py @@ -0,0 +1,313 @@ +from __future__ import annotations + +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 + +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)) + 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) + + # 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 + ) + + +@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) + laws = [ + ( + LinearElasticityModel( + parameters={"E": youngs_modulus, "nu": poissons_ratio}, + constraint=Constraint.UNIAXIAL_STRESS, + ), + np.array([0], dtype=np.int32), + ), + ( + LinearElasticityModel( + parameters={"E": factor * 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() + + # 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] - 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 + ) + + +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() + 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 + 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 + ) diff --git a/pyproject.toml b/pyproject.toml index 67d0ed9..9dde59e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -74,6 +74,7 @@ filterwarnings = [ log_cli_level = "INFO" testpaths = [ "tests", + "examples", ] @@ -102,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 a41c97c..2dc7c6b 100644 --- a/src/fenics_constitutive/__init__.py +++ b/src/fenics_constitutive/__init__.py @@ -9,6 +9,8 @@ from ._version import version as __version__ from .interfaces 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 new file mode 100644 index 0000000..e84e06e --- /dev/null +++ b/src/fenics_constitutive/conversion.py @@ -0,0 +1,81 @@ +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. + + 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 _: + msg = "Only full constraint implemented" + raise NotImplementedError(msg) + return strain diff --git a/src/fenics_constitutive/interfaces.py b/src/fenics_constitutive/interfaces.py index f490661..7a2c550 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 @@ -79,9 +81,10 @@ 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. + """ + Evaluate the constitutive model and overwrite the stress, tangent and history. Args: del_t : The time increment. @@ -90,7 +93,6 @@ def evaluate( tangent : The tangent. history : The history variable(s). """ - pass @abstractmethod def update(self) -> None: @@ -100,7 +102,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 +111,6 @@ def constraint(self) -> Constraint: Returns: The constraint. """ - pass @property def stress_strain_dim(self) -> int: @@ -133,7 +133,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 @@ -144,4 +144,3 @@ def history_dim(self) -> int | dict[str, int | tuple[int, int]]: Returns: The dimension of the history variable(s). """ - pass diff --git a/src/fenics_constitutive/maps.py b/src/fenics_constitutive/maps.py new file mode 100644 index 0000000..b9420a8 --- /dev/null +++ b/src/fenics_constitutive/maps.py @@ -0,0 +1,122 @@ +from __future__ import annotations + +from dataclasses import dataclass + +import dolfinx as df +import numpy as np + +__all__ = ["SubSpaceMap", "build_subspace_map"] + + +@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 + ), "Parent mesh does not match" + assert sub.ufl_shape == parent.ufl_shape, "Shapes do not match" + + 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: + """ + 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 + ), "Parent mesh does not match" + assert sub.ufl_shape == parent.ufl_shape, "Shapes do not match" + + size = sub.ufl_element().value_size() + + parent_array = parent.x.array.reshape(-1, size) + sub_array = sub.x.array.reshape(-1, size) + sub_array[self.child] = parent_array[self.parent] + sub.x.scatter_forward() + + +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] +): + """ + 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() + 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 diff --git a/src/fenics_constitutive/solver.py b/src/fenics_constitutive/solver.py new file mode 100644 index 0000000..4e7b263 --- /dev/null +++ b/src/fenics_constitutive/solver.py @@ -0,0 +1,331 @@ +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 + + +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. + + """ + match law.history_dim: + case int(): + Qh = ufl.VectorElement( + "Quadrature", + mesh.ufl_cell(), + q_degree, + quad_scheme="default", + dim=law.history_dim, + ) + history_space = df.fem.FunctionSpace(mesh, 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", + mesh.ufl_cell(), + q_degree, + quad_scheme="default", + dim=value, + ) + elif isinstance(value, tuple): + Qh = ufl.TensorElement( + "Quadrature", + mesh.ufl_cell(), + q_degree, + quad_scheme="default", + shape=value, + ) + history_space = df.fem.FunctionSpace(mesh, Qh) + history[key] = df.fem.Function(history_space) + return 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, + u: df.fem.Function, + bcs: list[df.fem.DirichletBCMetaClass], + q_degree: int = 1, + form_compiler_options: dict | None = None, + jit_options: dict | None = None, + ): + mesh = u.function_space.mesh + 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( + 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(), + 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()), + ) + 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_grad_u = [] + 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"): + 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._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) + 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)) + + # 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) + 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(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_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._u = u + self._u0 = u.copy() + 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) + + self.del_grad_u_expr = df.fem.Expression( + ufl.nabla_grad(self._u - self._u0), self.q_points + ) + + @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 + 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 + + 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, but here + we use it to update the stress, tangent and history. + + Args: + x: The vector containing the latest solution + + """ + super().form(x) + + assert ( + 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_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]) + 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 + if law.history_dim is not None + else None, + ) + + 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"): + 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 + 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[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, + ) + + self.stress_1.x.scatter_forward() + self.tangent.x.scatter_forward() + + def update(self) -> None: + """ + Update the current displacement, stress and history. + """ + 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/src/fenics_constitutive/stress_strain.py b/src/fenics_constitutive/stress_strain.py index c168607..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,9 @@ def ufl_mandel_strain( Returns: 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) @@ -132,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 a6aa8de..94b54a0 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(): @@ -41,26 +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) - 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) @@ -73,13 +83,9 @@ def test_ufl_strain_equals_array_conversion(constraint: Constraint): 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 new file mode 100644 index 0000000..0ab0c4e --- /dev/null +++ b/tests/test_maps.py @@ -0,0 +1,117 @@ +from __future__ import annotations + +import dolfinx as df +import numpy as np +import ufl +from mpi4py import MPI + +from fenics_constitutive import build_subspace_map + + +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) + 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_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) + + for _ in range(10): + cell_sample = np.random.permutation(cells)[: num_cells // 2] + + 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) + 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() + + 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] + + 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 = 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 = 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 = 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])