diff --git a/README.md b/README.md index a86152dd..1c0c1a51 100644 --- a/README.md +++ b/README.md @@ -13,7 +13,7 @@ Install Python, fire up a terminal and run pip install felupe[all] ``` -where `[all]` installs all optional dependencies. By default, FElupe does not require `numba` and `sparse`. In order to make use of all features of FElupe, it is suggested to install all optional dependencies. +where `[all]` installs all optional dependencies. By default, FElupe only depends on `numpy` and `scipy`. In order to make use of all features of FElupe, it is suggested to install all optional dependencies. # Hello, FElupe! A quarter model of a solid cube with hyperelastic material behavior is subjected to a uniaxial elongation applied at a clamped end-face. This involves the creation of a mesh, a region and a displacement field. Furthermore, the boundary conditions are created by a template for a uniaxial loadcase. The material behavior is defined through a FElupe-built-in Neo-Hookean material formulation. Inside a Newton-Rhapson procedure, the internal force vector and the tangent stiffness matrix are generated by assembling both linear and bilinear forms of static equilibrium. Finally, the solution of the incremental displacements is calculated und updated until convergence is reached. For more details beside this high-level code snippet, please have a look at the [documentation](https://felupe.readthedocs.io/en/latest/?badge=latest). @@ -49,7 +49,8 @@ All notable changes to this project will be documented in this file. The format ## [Unreleased] ## Added -Add optional parallel (threaded) basis evaluation and add `Form(v, u, parallel=True)`. +- Add optional parallel (threaded) basis evaluation and add `Form(v, u, parallel=True)`. +- Add `mechanics` submodule with `SolidBody` and `SolidBodyPressure`. ## [3.0.0] - 2022-04-28 diff --git a/docs/felupe/global.rst b/docs/felupe/global.rst index e829bf79..21618f4f 100644 --- a/docs/felupe/global.rst +++ b/docs/felupe/global.rst @@ -14,4 +14,5 @@ felupe global/dof global/constitution global/assembly - global/tools \ No newline at end of file + global/tools + global/mechanics \ No newline at end of file diff --git a/docs/felupe/global/mechanics.rst b/docs/felupe/global/mechanics.rst new file mode 100644 index 00000000..8ce17172 --- /dev/null +++ b/docs/felupe/global/mechanics.rst @@ -0,0 +1,12 @@ +Element +~~~~~~~ + +.. autoclass:: felupe.SolidBody + :members: + :undoc-members: + :show-inheritance: + +.. autoclass:: felupe.SolidBodyPressure + :members: + :undoc-members: + :show-inheritance: \ No newline at end of file diff --git a/docs/howto.rst b/docs/howto.rst index b733421e..b2d51093 100644 --- a/docs/howto.rst +++ b/docs/howto.rst @@ -12,3 +12,4 @@ How-To Guides howto/solvers howto/composite howto/forms + howto/solid \ No newline at end of file diff --git a/docs/howto/solid.rst b/docs/howto/solid.rst new file mode 100644 index 00000000..9a1a2e04 --- /dev/null +++ b/docs/howto/solid.rst @@ -0,0 +1,62 @@ +Solid Body (Mechanics) +~~~~~~~~~~~~~~~~~~~~~~ + +The generation of internal force vectors or stiffness matrices of solid bodies are provided as assembly-methods of a :class:`felupe.SolidBody`. The correct integral form is chosen based on the :class:`felupe.Field` (default, axisymmetric or mixed). + +.. code-block:: python + + import felupe as fe + + neohooke = fe.NeoHooke(mu=1.0, bulk=5000.0) + mesh = fe.Cube(n=6) + region = fe.RegionHexahedron(mesh) + displacement = fe.Field(region, dim=3) + + body = fe.SolidBody(umat=neohooke, field=displacement) + internal_force = body.assemble.vector(displacement, parallel=False, jit=False) + stiffness_matrix = body.assemble.matrix(displacement, parallel=False, jit=False) + + +During assembly, several results are stored, e.g. the gradient of the strain energy density function per unit undeformed volume (first Piola-Kirchhoff stress tensor). Other results are the deformation gradient or the fourth-order elasticity tensor associated to the first Piola-Kirchhoff stress tensor. + +.. code-block:: python + + F = body.results.kinematics[0] + P = body.results.stress + A = body.results.elasticity + + +The Cauchy stress tensor, as well as the gradient and the hessian of the strain energy density function per unit undeformed volume are obtained by evaluate-methods of the solid body. + +.. code-block:: python + + P = body.evaluate.gradient(displacement) + A = body.evaluate.hessian(displacement) + s = body.evaluate.cauchy_stress(displacement) + + +Pressure Boundary on Solid Body (Mechanics) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The generation of internal force vectors or stiffness matrices of pressure boundaries of solid bodies are provided as assembly-methods of a :class:`felupe.SolidBodyPressure`. The correct integral form is chosen based on the :class:`felupe.Field` (default or axisymmetric). If the internal field is a mixed field, the assembled vectors and matrices from the pressure contribution have to be resized to the dimensions of the internal force vector and the stiffness matrix. + +.. code-block:: python + + region_pressure = fe.RegionHexahedronBoundary( + mesh=mesh, + only_surface=True, # select only faces on the outline + mask=lambda x: x==0, # select a subset of faces on the surface + ) + + displacement_boundary = fe.Field(region_pressure) + displacement_boundary.values = displacement.values # link field values + + body_pressure = fe.SolidBodyPressure(field=displacement) + + internal_force_pressure = body.assemble.vector( + field=displacement, parallel=False, jit=False, resize=internal_force + ) + + stiffness_matrix_pressure = body.assemble.matrix( + field=displacement, parallel=False, jit=False, resize=stiffness_matrix + ) \ No newline at end of file diff --git a/docs/howto/solvers.rst b/docs/howto/solvers.rst index 0a772b56..055e7296 100644 --- a/docs/howto/solvers.rst +++ b/docs/howto/solvers.rst @@ -88,36 +88,15 @@ Solvers from external packages: # ... def solver(A, b): - return PyPardisoSolver(mtype=6).solve(triu(A).tocsr(), b).squeeze() + # mtype = 1: real and structurally symmetric, supernode pivoting + # mtype = 2: real and symmetric positive definite + # mtype =-2: real and symmetric indefinite, + # diagonal or Bunch-Kaufman pivoting + # mtype = 6: complex and symmetric + return PyPardisoSolver(mtype=2).solve(triu(A).tocsr(), b).squeeze() system = fe.solve.partition(field, K, dof1, dof0) fe.solve.solve(*system, solver=solver) -.. tab:: Krylov (iterative) - - Ensure to have `Krylov `_ installed. - - .. code-block:: bash - - pip install krylov - - ``minres`` may be replaced by another iterative method. - - .. code-block:: python - - import felupe as fe - import krylov - - # ... - - def solver(A, b): - "Wrapper function for Krylov-solvers." - - return krylov.minres(A, b)[0] - - system = fe.solve.partition(field, K, dof1, dof0) - fe.solve.solve(*system, solver=solver) - - \ No newline at end of file diff --git a/docs/index.rst b/docs/index.rst index bc85dedc..b7f08e60 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -57,7 +57,7 @@ Another key feature is the easy and straightforward definition of mixed field fo Installation ------------ -Install Python, fire up a terminal and run ``pip install felupe[all]``, where ``[all]`` installs all optional dependencies. By default, FElupe only depends on ``numpy``, ``scipy`` and ``meshio``. However, ``h5py``, ``numba``, ``einsumt`` and ``sparse`` are highly recommended. In order to make use of all features of FElupe, it is suggested to install all optional dependencies. For constitutive material definitions using Automatic Differentation please also install `matADi `_. +Install Python, fire up a terminal and run ``pip install felupe[all]``, where ``[all]`` installs all optional dependencies. By default, FElupe only depends on ``numpy`` and ``scipy``. However, ``meshio``, ``h5py``, ``numba``, ``einsumt`` and ``sparse`` are highly recommended. In order to make use of all features of FElupe, it is suggested to install all optional dependencies. For constitutive material definitions using Automatic Differentation consider also installing `matADi `_. .. code-block:: shell diff --git a/felupe/__about__.py b/felupe/__about__.py index 528787cf..f5f41e56 100644 --- a/felupe/__about__.py +++ b/felupe/__about__.py @@ -1 +1 @@ -__version__ = "3.0.0" +__version__ = "3.1.0" diff --git a/felupe/__init__.py b/felupe/__init__.py index 57791d4b..5c5f382b 100644 --- a/felupe/__init__.py +++ b/felupe/__init__.py @@ -8,6 +8,7 @@ from . import constitution from . import solve from . import region +from . import mechanics from .region import ( Region, @@ -107,6 +108,10 @@ topoints, project, ) +from .mechanics import ( + SolidBody, + SolidBodyPressure, +) __all__ = [ "__version__", diff --git a/felupe/_assembly/_axi.py b/felupe/_assembly/_axi.py index 2187efff..acf6a05d 100644 --- a/felupe/_assembly/_axi.py +++ b/felupe/_assembly/_axi.py @@ -46,10 +46,10 @@ def __init__(self, fun, v, dV, u=None, grad_v=True, grad_u=True): if grad_v: fun_2d = fun[:-1, :-1] - fun_zz = fun[(-1,), (-1,)] + fun_zz = fun[(-1,), (-1,)] / R else: fun_2d = fun[:-1] - fun_zz = fun[-1].reshape(1, *fun[-1].shape) + fun_zz = fun[-1].reshape(1, *fun[-1].shape) / R form_a = IntegralForm(fun_2d, v, self.dV, grad_v=grad_v) form_b = IntegralForm(fun_zz, v.scalar, self.dV) @@ -142,6 +142,10 @@ def integrate(self, parallel=False, jit=False): val = values[0] if self.mode == 30: + if len(values[0].shape) > 4: + values[0] = values[0][:, :, 0, 0] + if len(values[1].shape) > 4: + values[1] = values[1][:, :, 0, 0] a, b, e = values[1].shape values[1] = values[1].reshape(a, 1, b, e) diff --git a/felupe/constitution/_kinematics.py b/felupe/constitution/_kinematics.py index 97f69db8..53fda599 100644 --- a/felupe/constitution/_kinematics.py +++ b/felupe/constitution/_kinematics.py @@ -77,7 +77,7 @@ def function(self, F): """ return F - def gradient(self, F): + def gradient(self, F, parallel=None): """Gradient of line change. Arguments @@ -91,8 +91,11 @@ def gradient(self, F): Gradient of line change """ + if parallel is None: + parallel = self.parallel + Eye = identity(F) - return cdya_ik(Eye, Eye, parallel=self.parallel) + return cdya_ik(Eye, Eye, parallel=parallel) class AreaChange: @@ -113,7 +116,7 @@ class AreaChange: def __init__(self, parallel=False): self.parallel = parallel - def function(self, F, N=None): + def function(self, F, N=None, parallel=None): """Area change. Arguments @@ -132,12 +135,15 @@ def function(self, F, N=None): Fs = J * transpose(inv(F, J)) + if parallel is None: + parallel = self.parallel + if N is None: return Fs else: - return dot(Fs, N, parallel=self.parallel) + return dot(Fs, N, parallel=parallel) - def gradient(self, F, N=None): + def gradient(self, F, N=None, parallel=None): """Gradient of area change. Arguments @@ -155,13 +161,15 @@ def gradient(self, F, N=None): J = det(F) + if parallel is None: + parallel = self.parallel + dJdF = self.function(F) dFsdF = ( - dya(dJdF, dJdF, parallel=self.parallel) - - cdya_il(dJdF, dJdF, parallel=self.parallel) + dya(dJdF, dJdF, parallel=parallel) - cdya_il(dJdF, dJdF, parallel=parallel) ) / J - if self.parallel: + if parallel: einsum = einsumt else: einsum = np.einsum @@ -224,7 +232,7 @@ def gradient(self, F): J = self.function(F) return J * transpose(inv(F, J)) - def hessian(self, F): + def hessian(self, F, parallel=None): """Hessian of volume change. Arguments @@ -238,9 +246,11 @@ def hessian(self, F): Hessian of the determinant of the deformation gradient """ + if parallel is None: + parallel = self.parallel + J = self.function(F) dJdF = self.gradient(F) return ( - dya(dJdF, dJdF, parallel=self.parallel) - - cdya_il(dJdF, dJdF, parallel=self.parallel) + dya(dJdF, dJdF, parallel=parallel) - cdya_il(dJdF, dJdF, parallel=parallel) ) / J diff --git a/felupe/mechanics/__init__.py b/felupe/mechanics/__init__.py new file mode 100644 index 00000000..81657ed6 --- /dev/null +++ b/felupe/mechanics/__init__.py @@ -0,0 +1,2 @@ +from ._solidbody import SolidBody +from ._solidbody_pressure import SolidBodyPressure diff --git a/felupe/mechanics/_helpers.py b/felupe/mechanics/_helpers.py new file mode 100644 index 00000000..6c7f4339 --- /dev/null +++ b/felupe/mechanics/_helpers.py @@ -0,0 +1,61 @@ +# -*- coding: utf-8 -*- +""" + _______ _______ ___ __ __ _______ _______ +| || || | | | | || || | +| ___|| ___|| | | | | || _ || ___| +| |___ | |___ | | | |_| || |_| || |___ +| ___|| ___|| |___ | || ___|| ___| +| | | |___ | || || | | |___ +|___| |_______||_______||_______||___| |_______| + +This file is part of felupe. + +Felupe is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +Felupe is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with Felupe. If not, see . + +""" + + +class Assemble: + "A class with assembly methods of a SolidBody." + + def __init__(self, vector, matrix): + self.vector = vector + self.matrix = matrix + + +class Evaluate: + "A class with evaluate methods of a SolidBody." + + def __init__(self, gradient, hessian, cauchy_stress=None): + self.gradient = gradient + self.hessian = hessian + + if cauchy_stress is not None: + self.cauchy_stress = cauchy_stress + + +class Results: + "A class with intermediate results of a SolidBody." + + def __init__(self, stress=False, elasticity=False): + + self.force = None + self.stiffness = None + self.kinematics = None + + if stress: + self.stress = None + + if elasticity: + self.elasticity = None diff --git a/felupe/mechanics/_solidbody.py b/felupe/mechanics/_solidbody.py new file mode 100644 index 00000000..707eec51 --- /dev/null +++ b/felupe/mechanics/_solidbody.py @@ -0,0 +1,155 @@ +# -*- coding: utf-8 -*- +""" + _______ _______ ___ __ __ _______ _______ +| || || | | | | || || | +| ___|| ___|| | | | | || _ || ___| +| |___ | |___ | | | |_| || |_| || |___ +| ___|| ___|| |___ | || ___|| ___| +| | | |___ | || || | | |___ +|___| |_______||_______||_______||___| |_______| + +This file is part of felupe. + +Felupe is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +Felupe is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with Felupe. If not, see . + +""" + +from .._field import Field, FieldMixed, FieldsMixed, FieldAxisymmetric +from .._assembly import IntegralForm, IntegralFormMixed, IntegralFormAxisymmetric +from ..constitution import AreaChange +from ..math import inv, dot +from ._helpers import Assemble, Evaluate, Results + + +class SolidBody: + "A SolidBody with methods for the assembly of sparse vectors/matrices." + + def __init__(self, umat, field): + + self.umat = umat + self.field = field + + if isinstance(field, FieldMixed): + self._dV = self.field[0].region.dV + else: + self._dV = self.field.region.dV + + self.results = Results(stress=True, elasticity=True) + self.results.kinematics = self._extract(self.field) + + self.assemble = Assemble(vector=self._vector, matrix=self._matrix) + + self.evaluate = Evaluate( + gradient=self._gradient, + hessian=self._hessian, + cauchy_stress=self._cauchy_stress, + ) + + self._area_change = AreaChange() + + self._form = { + Field: IntegralForm, + FieldMixed: IntegralFormMixed, + FieldsMixed: IntegralFormMixed, + FieldAxisymmetric: IntegralFormAxisymmetric, + }[type(self.field)] + + self._kwargs = { + Field: dict(dV=self._dV, grad_v=True, grad_u=True), + FieldMixed: dict(dV=self._dV), + FieldsMixed: dict(dV=self._dV), + FieldAxisymmetric: dict(dV=self._dV, grad_v=True, grad_u=True), + }[type(self.field)] + + def _vector( + self, field=None, parallel=False, jit=False, items=None, args=(), kwargs={} + ): + + if field is not None: + self.field = field + + self.results.stress = self._gradient(field, args=args, kwargs=kwargs) + + self.results.force = self._form( + fun=self.results.stress[slice(items)], + v=self.field, + **self._kwargs, + ).assemble(parallel=parallel, jit=jit) + + return self.results.force + + def _matrix( + self, field=None, parallel=False, jit=False, items=None, args=(), kwargs={} + ): + + if field is not None: + self.field = field + + self.results.elasticity = self._hessian(field, args=args, kwargs=kwargs) + + self.results.stiffness = self._form( + fun=self.results.elasticity[slice(items)], + v=self.field, + u=self.field, + **self._kwargs, + ).assemble(parallel=parallel, jit=jit) + + return self.results.stiffness + + def _extract(self, field): + + self.field = field + + self.results.kinematics = self.field.extract() + if isinstance(self.field, Field): + self.results.kinematics = (self.results.kinematics,) + + return self.results.kinematics + + def _gradient(self, field=None, args=(), kwargs={}): + + if field is not None: + self.field = field + self.results.kinematics = self._extract(self.field) + + self.results.stress = self.umat.gradient( + *self.results.kinematics, *args, **kwargs + ) + + return self.results.stress + + def _hessian(self, field=None, args=(), kwargs={}): + + if field is not None: + self.field = field + self.results.kinematics = self._extract(self.field) + + self.results.elasticity = self.umat.hessian( + *self.results.kinematics, *args, **kwargs + ) + + return self.results.elasticity + + def _cauchy_stress(self, field=None): + + self._gradient(field) + + if len(self.results.kinematics) > 1: + P = self.results.stress[0] + else: + P = self.results.stress + + JiFT = self._area_change.function(self.results.kinematics[0]) + + return dot(P, inv(JiFT)) diff --git a/felupe/mechanics/_solidbody_pressure.py b/felupe/mechanics/_solidbody_pressure.py new file mode 100644 index 00000000..2af111d7 --- /dev/null +++ b/felupe/mechanics/_solidbody_pressure.py @@ -0,0 +1,138 @@ +# -*- coding: utf-8 -*- +""" + _______ _______ ___ __ __ _______ _______ +| || || | | | | || || | +| ___|| ___|| | | | | || _ || ___| +| |___ | |___ | | | |_| || |_| || |___ +| ___|| ___|| |___ | || ___|| ___| +| | | |___ | || || | | |___ +|___| |_______||_______||_______||___| |_______| + +This file is part of felupe. + +Felupe is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +Felupe is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with Felupe. If not, see . + +""" + +from .._field import Field, FieldMixed, FieldsMixed, FieldAxisymmetric +from .._assembly import IntegralForm, IntegralFormMixed, IntegralFormAxisymmetric +from ..constitution import AreaChange +from ._helpers import Assemble, Results + + +class SolidBodyPressure: + "A hydrostatic pressure boundary on a SolidBody." + + def __init__(self, field): + + self.field = field + + self._dV = self.field.region.dV + self._normals = self.field.region.normals + + self.results = Results() + self.results.kinematics = self._extract(self.field) + self.assemble = Assemble(vector=self._vector, matrix=self._matrix) + + self._form = { + Field: IntegralForm, + FieldMixed: IntegralFormMixed, + FieldsMixed: IntegralFormMixed, + FieldAxisymmetric: IntegralFormAxisymmetric, + }[type(self.field)] + + self._kwargs = { + Field: dict(dV=self._dV, grad_v=True, grad_u=True), + FieldMixed: dict(dV=self._dV), + FieldsMixed: dict(dV=self._dV), + FieldAxisymmetric: dict(dV=self._dV, grad_v=True, grad_u=True), + }[type(self.field)] + + self._IntForm = { + Field: IntegralForm, + FieldAxisymmetric: IntegralFormAxisymmetric, + }[type(self.field)] + + self._area_change = AreaChange() + + def _extract(self, field): + + self.field = field + self.results.kinematics = (self.field.extract(),) + + return self.results.kinematics + + def _vector(self, field=None, pressure=1, parallel=False, jit=False, resize=None): + + if field is not None: + self.field = field + self.results.kinematics = self._extract(field) + + self.results.pressure = pressure + + fun = pressure * self._area_change.function( + *self.results.kinematics, + self._normals, + parallel=parallel, + ) + + self.results.force = self._IntForm( + fun=fun, v=self.field, dV=self._dV, grad_v=False + ).assemble(parallel=parallel, jit=jit) + + if resize is not None: + self.results.force.resize(*resize.shape) + + return self.results.force + + def _matrix(self, field=None, pressure=1, parallel=False, jit=False, resize=None): + + if field is not None: + self.field = field + self.results.kinematics = self._extract(field) + + self.results.pressure = pressure + + fun = pressure * self._area_change.gradient( + *self.results.kinematics, + self._normals, + parallel=parallel, + ) + self.results.stiffness = self._IntForm( + fun=fun, + v=self.field, + u=self.field, + dV=self._dV, + grad_v=False, + grad_u=True, + ).assemble(parallel=parallel, jit=jit) + + if resize is not None: + self.results.stiffness.resize(*resize.shape) + + return self.results.stiffness + + def update(self, other_field, field=None): + + if field is not None: + self.field = field + + if isinstance(other_field, FieldMixed) or isinstance(other_field, FieldsMixed): + self.field.values = other_field[0].values + else: + self.field.values = other_field.values + + self.results.kinematics = self._extract(self.field) + + return self.field diff --git a/felupe/mesh/_mesh.py b/felupe/mesh/_mesh.py index a1d1fba3..d9326531 100644 --- a/felupe/mesh/_mesh.py +++ b/felupe/mesh/_mesh.py @@ -117,10 +117,7 @@ def as_meshio(self, **kwargs): """ - if self.cell_type is None: - raise TypeError("Cell type missing.") - else: - import meshio + import meshio cells = {self.cell_type: self.cells} return meshio.Mesh(self.points, cells, **kwargs) diff --git a/felupe/tools/_save.py b/felupe/tools/_save.py index 91193be7..0a927177 100644 --- a/felupe/tools/_save.py +++ b/felupe/tools/_save.py @@ -26,7 +26,6 @@ """ import numpy as np -import meshio from ..math import dot, transpose, det, eigvalsh from . import topoints @@ -90,6 +89,8 @@ def save( point_data["MaxPrincipalShearCauchyStress"] = cauchyprinc[2] - cauchyprinc[0] + import meshio + mesh = meshio.Mesh( points=mesh.points, cells=[ diff --git a/pyproject.toml b/pyproject.toml index c6078a48..e16166e7 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "felupe" -version = "3.0.0" +version = "3.1.0" description = "Finite Element Analysis" readme = "README.md" requires-python = ">=3.6" @@ -40,11 +40,11 @@ classifiers = [ dependencies = [ "numpy", "scipy", - "meshio", ] [project.optional-dependencies] all = [ + "meshio", "h5py", "numba", "sparse", diff --git a/setup.cfg b/setup.cfg index ebdd5a85..5a983652 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,6 +1,6 @@ [metadata] name = felupe -version = 3.0.0 +version = 3.1.0 author = Andreas Dutzler author_email = a.dutzler@gmail.com description = Finite Element Analysis @@ -33,11 +33,11 @@ packages = find: install_requires = numpy scipy - meshio python_requires = >=3.6 [options.extras_require] all = + meshio h5py numba sparse diff --git a/tests/test_constitution.py b/tests/test_constitution.py index 183c53d7..9b82777c 100644 --- a/tests/test_constitution.py +++ b/tests/test_constitution.py @@ -47,7 +47,7 @@ def pre_mixed(sym, add_identity): v = fe.Field(r, dim=1) z = fe.Field(r, dim=1, values=1) w = fe.FieldMixed((u, v, z)) - return r, w.extract(grad=True, sym=sym, add_identity=add_identity) + return r, w.extract(grad=True, sym=sym, add_identity=add_identity), m def test_nh(): @@ -207,9 +207,12 @@ def test_kinematics(): xf = lc.function(F) xg = lc.gradient(F) + xg = lc.gradient(F, parallel=parallel) Yf = ac.function(F, N) + Yf = ac.function(F, N, parallel=parallel) Yg = ac.gradient(F, N) + Yg = ac.gradient(F, N, parallel=parallel) yf = ac.function(F) yg = ac.gradient(F) @@ -217,6 +220,7 @@ def test_kinematics(): zf = vc.function(F) zg = vc.gradient(F) zh = vc.hessian(F) + zh = vc.hessian(F, parallel=parallel) assert np.allclose(xf, F) @@ -261,6 +265,13 @@ def hessian(self, x, threads=1): if len(x) == 1: return [self.material.hessian(*x)] else: + hess = self.material.hessian(*x) + for a in [1, 2]: + hess[a].reshape(*hess[a].shape[:-2], 1, 1, *hess[a].shape[-2:]) + for b in [3, 4, 5]: + hess[b].reshape( + *hess[b].shape[:-2], 1, 1, 1, *hess[b].shape[-2:] + ) return self.material.hessian(*x) umat = fe.MatadiMaterial(AsMatadi(nh)) @@ -273,7 +284,7 @@ def hessian(self, x, threads=1): assert P.shape == (3, 3, *F.shape[-2:]) assert A.shape == (3, 3, 3, 3, *F.shape[-2:]) - r, FpJ = pre_mixed(sym=False, add_identity=True) + r, FpJ, m = pre_mixed(sym=False, add_identity=True) umat = fe.MatadiMaterial( AsMatadi(fe.ThreeFieldVariation(nh, parallel=parallel)) @@ -285,6 +296,15 @@ def hessian(self, x, threads=1): assert P[0].shape == (3, 3, *FpJ[0].shape[-2:]) assert A[0].shape == (3, 3, 3, 3, *FpJ[0].shape[-2:]) + m = fe.Rectangle(n=3) + r = fe.RegionQuad(m) + v = fe.FieldsMixed(r, n=3, axisymmetric=True) + FpJ = v.extract() + A = umat.hessian(*FpJ) + K = fe.IntegralFormMixed(A, v, r.dV, v).assemble() + + assert K.shape == (26, 26) + if __name__ == "__main__": test_nh() diff --git a/tests/test_mechanics.py b/tests/test_mechanics.py new file mode 100644 index 00000000..1d9b2c62 --- /dev/null +++ b/tests/test_mechanics.py @@ -0,0 +1,289 @@ +# -*- coding: utf-8 -*- +""" + _______ _______ ___ __ __ _______ _______ +| || || | | | | || || | +| ___|| ___|| | | | | || _ || ___| +| |___ | |___ | | | |_| || |_| || |___ +| ___|| ___|| |___ | || ___|| ___| +| | | |___ | || || | | |___ +|___| |_______||_______||_______||___| |_______| + +This file is part of felupe. + +Felupe is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +Felupe is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with Felupe. If not, see . + +""" + +import pytest +import felupe as fe +import numpy as np + + +def test_simple(): + + umat = fe.LinearElastic(E=1, nu=0.3) + + m = fe.Cube(n=3) + r = fe.RegionHexahedron(m) + u = fe.Field(r, dim=3) + + b = fe.SolidBody(umat, u) + r = b.assemble.vector() + + K = b.assemble.matrix() + r = b.assemble.vector(u) + F = b.results.kinematics[0] + P = b.results.stress + s = b.evaluate.cauchy_stress() + C = b.results.elasticity + + assert K.shape == (81, 81) + assert r.shape == (81, 1) + assert F.shape == (3, 3, 8, 8) + assert P.shape == (3, 3, 8, 8) + assert s.shape == (3, 3, 8, 8) + assert C.shape == (3, 3, 3, 3, 8, 8) + + +def test_pressure(): + + umat = fe.LinearElastic(E=1, nu=0.3) + + m = fe.Cube(n=3) + h = fe.RegionHexahedron(m) + u = fe.Field(h, dim=3) + + u.values = np.random.rand(*u.values.shape) / 10 + + s = fe.RegionHexahedronBoundary(m) + v = fe.Field(s, dim=3) + + b = fe.SolidBody(umat, u) + c = fe.SolidBodyPressure(v) + + r = b.assemble.vector() + K = b.assemble.matrix() + r = b.assemble.vector(u) + F = b.results.kinematics[0] + s = b.results.stress + C = b.results.elasticity + + assert K.shape == (81, 81) + assert r.shape == (81, 1) + assert F.shape == (3, 3, 8, 8) + assert s.shape == (3, 3, 8, 8) + assert C.shape == (3, 3, 3, 3, 8, 8) + + r = c.assemble.vector() + K = c.assemble.matrix() + K = c.assemble.matrix(v, resize=b.assemble.matrix()) + r = c.assemble.vector(v) + r = c.assemble.vector(v, resize=b.assemble.vector()) + F = c.results.kinematics[0] + + assert K.shape == (81, 81) + assert r.shape == (81, 1) + assert F.shape == (3, 3, 4, 24) + + c.update(u, v) + assert np.allclose(u.values, v.values) + + w = fe.FieldsMixed(h) + w[0].values = np.random.rand(*w[0].values.shape) / 10 + + c.update(w, v) + assert np.allclose(w[0].values, v.values) + + +def pre(dim): + + umat = fe.NeoHooke(mu=1, bulk=2) + + m = fe.Cube(n=3) + r = fe.RegionHexahedron(m) + u = fe.Field(r, dim=dim) + + u.values = np.random.rand(*u.values.shape) / 10 + + return umat, u + + +def pre_axi(): + + umat = fe.NeoHooke(mu=1, bulk=2) + + m = fe.Rectangle(n=3) + r = fe.RegionQuad(m) + u = fe.FieldAxisymmetric(r) + + u.values = np.random.rand(*u.values.shape) / 10 + + return umat, u + + +def pre_mixed(dim): + + umat = fe.ThreeFieldVariation(fe.NeoHooke(mu=1, bulk=2)) + + m = fe.Cube(n=3) + r = fe.RegionHexahedron(m) + u = fe.FieldsMixed(r, n=3) + + u[0].values = np.random.rand(*u[0].values.shape) / 10 + u[1].values = np.random.rand(*u[1].values.shape) / 10 + u[2].values = np.random.rand(*u[2].values.shape) / 10 + + return umat, u + + +def test_solidbody(): + + umat, u = pre(dim=3) + b = fe.SolidBody(umat=umat, field=u) + + for parallel in [False, True]: + for jit in [False, True]: + + kwargs = {"parallel": parallel, "jit": jit} + + r1 = b.assemble.vector(u, **kwargs) + assert r1.shape == (81, 1) + + r1b = b.results.force + assert np.allclose(r1.toarray(), r1b.toarray()) + + r2 = b.assemble.vector(**kwargs) + assert np.allclose(r1.toarray(), r2.toarray()) + + K1 = b.assemble.matrix(u, **kwargs) + assert K1.shape == (81, 81) + + K1b = b.results.stiffness + assert np.allclose(K1.toarray(), K1b.toarray()) + + K2 = b.assemble.matrix(**kwargs) + assert np.allclose(K1.toarray(), K2.toarray()) + + P1 = b.results.stress + P2 = b.evaluate.gradient() + P2 = b.evaluate.gradient(u) + assert np.allclose(P1, P2) + + A1 = b.results.elasticity + A2 = b.evaluate.hessian() + A2 = b.evaluate.hessian(u) + assert np.allclose(A1, A2) + + F1 = b.results.kinematics + F2 = b._extract(u) + assert np.allclose(F1, F2) + + s1 = b.evaluate.cauchy_stress() + s2 = b.evaluate.cauchy_stress(u) + assert np.allclose(s1, s2) + + +def test_solidbody_axi(): + + umat, u = pre_axi() + b = fe.SolidBody(umat=umat, field=u) + + for parallel in [False, True]: + for jit in [False, True]: + + kwargs = {"parallel": parallel, "jit": jit} + + r1 = b.assemble.vector(u, **kwargs) + assert r1.shape == (18, 1) + + r2 = b.assemble.vector(**kwargs) + assert np.allclose(r1.toarray(), r2.toarray()) + + K1 = b.assemble.matrix(u, **kwargs) + assert K1.shape == (18, 18) + + K2 = b.assemble.matrix(**kwargs) + assert np.allclose(K1.toarray(), K2.toarray()) + + P1 = b.results.stress + P2 = b.evaluate.gradient() + P2 = b.evaluate.gradient(u) + assert np.allclose(P1, P2) + + A1 = b.results.elasticity + A2 = b.evaluate.hessian() + A2 = b.evaluate.hessian(u) + assert np.allclose(A1, A2) + + F1 = b.results.kinematics + F2 = b._extract(u) + assert np.allclose(F1, F2) + + s1 = b.evaluate.cauchy_stress() + s2 = b.evaluate.cauchy_stress(u) + assert np.allclose(s1, s2) + + +def test_solidbody_mixed(): + + umat, u = pre_mixed(dim=3) + b = fe.SolidBody(umat=umat, field=u) + + for parallel in [False, True]: + for jit in [False, True]: + + kwargs = {"parallel": parallel, "jit": jit} + + r1 = b.assemble.vector(u, **kwargs) + r1 = b.assemble.vector(u, items=3, **kwargs) + assert r1.shape == (97, 1) + + r2 = b.assemble.vector(**kwargs) + assert np.allclose(r1.toarray(), r2.toarray()) + + K1 = b.assemble.matrix(u, **kwargs) + K1 = b.assemble.matrix(u, items=6, **kwargs) + assert K1.shape == (97, 97) + + K2 = b.assemble.matrix(**kwargs) + assert np.allclose(K1.toarray(), K2.toarray()) + + P1 = b.results.stress + P2 = b.evaluate.gradient() + P2 = b.evaluate.gradient(u) + for p1, p2 in zip(P1, P2): + assert np.allclose(p1, p2) + + A1 = b.results.elasticity + A2 = b.evaluate.hessian() + A2 = b.evaluate.hessian(u) + for a1, a2 in zip(A1, A2): + assert np.allclose(a1, a2) + + F1 = b.results.kinematics + F2 = b._extract(u) + for f1, f2 in zip(F1, F2): + assert np.allclose(f1, f2) + + s1 = b.evaluate.cauchy_stress() + s2 = b.evaluate.cauchy_stress(u) + assert np.allclose(s1, s2) + + +if __name__ == "__main__": + test_simple() + test_solidbody() + test_solidbody_axi() + test_solidbody_mixed() + test_pressure() diff --git a/tests/test_mesh.py b/tests/test_mesh.py index 4c9df968..695fd7e5 100644 --- a/tests/test_mesh.py +++ b/tests/test_mesh.py @@ -142,7 +142,7 @@ def test_meshes(): m.save() m.cell_type = None - with pytest.raises(TypeError): + with pytest.raises(Exception): m.save()