diff --git a/doc/dev-docs/source/ConstitutiveLaws.rst b/doc/dev-docs/source/ConstitutiveLaws.rst index 309b52b..9b6bc6a 100644 --- a/doc/dev-docs/source/ConstitutiveLaws.rst +++ b/doc/dev-docs/source/ConstitutiveLaws.rst @@ -1,15 +1,128 @@ .. toctree:: :maxdepth: 2 :caption: Contents: .. _constitutive_laws: Constitutive Laws ~~~~~~~~~~~~~~~~~ .. toctree:: :maxdepth: 2 MaterialLinearElasticGeneric + + +Testing Constitutive Laws +~~~~~~~~~~~~~~~~~~~~~~~~~ + +When writing new constitutive laws, the ability to evaluate the stress-strain +behaviour and the tangent moduli is convenient, but *µ*\Spectre's material model +makes it cumbersome to isolate and execute the :cpp:func:`evaluate_stress` and +:cpp:func:`evaluate_stress_tangent` methods than any daughter class of +:cpp:class:`MaterialMuSpectre` must implement +(e.g., :cpp:func:`evaluate_stress() +`). As a helper object, +*µ*\Spectre offers the class +:cpp:class:`MaterialEvaluator` to facilitate +precisely this: + +A :cpp:class:`MaterialEvaluator` object can be +constructed with a shared pointer to a +:cpp:class:`MaterialBase` and exposes functions to +evaluate just the stress, both the stress and tangent moduli, or a numerical +approximation to the tangent moduli. For materials with internal history +variables, :cpp:class:`MaterialEvaluator` also +exposes the +:cpp:func:`MaterialBase::save_history_variables()` +method. As a convenience function, all daughter classes of +:cpp:class:`MaterialMuSpectre` have the static +factory function :cpp:func:`make_evaluator() +` to create a material and its +evaluator at once. See the :ref:`reference` for the full class description. + +Python Usage Example +```````````````````` + +.. code-block:: python + + import numpy as np + from muSpectre import material + from muSpectre import Formulation + + # MaterialLinearElastic1 is standard linear elasticity while + # MaterialLinearElastic2 has a per pixel eigenstrain which needs to be set + LinMat1, LinMat2 = (material.MaterialLinearElastic1_2d, + material.MaterialLinearElastic2_2d) + + young, poisson = 210e9, .33 + + # the factory returns a material and it's corresponding evaluator + material1, evaluator1 = LinMat1.make_evaluator(young, poisson) + + # the material is empty (i.e., does not have any pixel/voxel), so a pixel + # needs to be added. The coordinates are irrelevant, there just needs to + # be one pixel. + material1.add_pixel([0,0]) + + # the stress and tangent can be evaluated for finite strain + F = np.array([[1., .01],[0, 1.0]]) + P, K = evaluator1.evaluate_stress_tangent(F, Formulation.finite_strain) + # or small strain + eps = .5 * ((F-np.eye(2)) + (F-np.eye(2)).T) + sigma, C = evaluator1.evaluate_stress_tangent(eps, Formulation.small_strain) + + # and the tangent can be checked against a numerical approximation + Delta_x = 1e-6 + num_C = evaluator1.estimate_tangent(eps, Formulation.small_strain, Delta_x) + + + # Materials with per-pixel data behave similarly: the factory returns a + # material and it's corresponding evaluator like before + material2, evaluator2 = LinMat2.make_evaluator(young, poisson) + + # when adding the pixel, we now need to specify also the per-pixel data: + eigenstrain = np.array([[.01, .002], [.002, 0.]]) + material2.add_pixel([0,0], eigenstrain) + + + +C++ Usage Example +````````````````` + +.. code-block:: c++ + + #include"materials/material_linear_elastic2.hh" + #include "materials/material_evaluator.hh" + #include "common/T4_map_proxy.hh" + + #include "Eigen/Dense" + + using Mat_t = MaterialLinearElastic2; + + constexpr Real Young{210e9}; + constexpr Real Poisson{.33}; + + auto mat_eval{Mat_t::make_evaluator(Young, Poisson)}; + auto & mat{*std::get<0>(mat_eval)}; + auto & evaluator{std::get<1>(mat_eval)}; + + + using T2_t = Eigen::Matrix; + using T4_t = T4Mat; + const T2_t F{(T2_t::Random() - (T2_t::Ones() * .5)) * 1e-4 + + T2_t::Identity()}; + + T2_t eigen_strain{[](auto x) { + return 1e-4 * (x + x.transpose()); + }(T2_t::Random() - T2_t::Ones() * .5)}; + + mat.add_pixel({}, eigen_strain); + + T2_t P{}; + T4_t K{}; + + std::tie(P, K) = + evaluator.evaluate_stress_tangent(F, Formulation::finite_strain); diff --git a/doc/dev-docs/source/Doxyfile b/doc/dev-docs/source/Doxyfile index 3b8c664..cd8fe05 100644 --- a/doc/dev-docs/source/Doxyfile +++ b/doc/dev-docs/source/Doxyfile @@ -1,28 +1,29 @@ PROJECT_NAME = muSpectre ## this is for read the docs, which builds in-place INPUT = ../../../src ## this overwrites the INPUT path for local builds @INCLUDE = input_def RECURSIVE = YES EXTRACT_ALL = YES SOURCE_BROWSER = YES GENERATE_HTML = YES GENERATE_LATEX = NO GENERATE_XML = YES ## this is for read the docs, which builds in-place XML_OUTPUT = doxygenxml ## this overwrites the XML_OUTPUT path for local builds @INCLUDE = xml_output_def XML_PROGRAMLISTING = YES ALIASES = "rst=\verbatim embed:rst" ALIASES += "endrst=\endverbatim" QUIET = YES WARNINGS = YES WARN_IF_UNDOCUMENTED = YES +SHOW_NAMESPACES = NO EXCLUDE_PATTERNS = */CMakeLists.txt \ No newline at end of file diff --git a/doc/dev-docs/source/GettingStarted.rst b/doc/dev-docs/source/GettingStarted.rst index 8e3007e..fb2eaa8 100644 --- a/doc/dev-docs/source/GettingStarted.rst +++ b/doc/dev-docs/source/GettingStarted.rst @@ -1,102 +1,136 @@ Getting Started ~~~~~~~~~~~~~~~ Obtaining *µ*\Spectre ********************* *µ*\Spectre is hosted on a git repository on `c4science`_. To clone it, run .. code-block:: sh $ git clone https://c4science.ch/source/muSpectre.git or if you prefer identifying yourself using a public ssh-key, run .. code-block:: bash $ git clone ssh://git@c4science.ch/source/muSpectre.git -The latter option requires you to have a user account on c4science (`create `_). +The latter option requires you to have a user account on c4science (`create +`_). .. _c4science: https://c4science.ch Building *µ*\Spectre ******************** -You can compile *µ*\Spectre using `CMake `_ (3.1.0 or higher). The current (and possibly incomplete list of) dependencies are +You can compile *µ*\Spectre using `CMake `_ (3.1.0 or +higher). The current (and possibly incomplete list of) dependencies are - `CMake `_ (3.1.0 or higher) - `Boost unit test framework `_ - `FFTW `_ - `git `_ - `Python3 `_ including the header files. - `numpy `_ and `scipy `_. Recommended: -- `Sphinx `_ and `Breathe `_ (necessary if you want to build the documentation (turned off by default) -- `Eigen `_ (3.3.0 or higher). If you do not install this, it will be downloaded automatically at configuration time, so this is not strictly necessary. The download can be slow, though, so we recommend installing it on your system. +- `Sphinx `_ and `Breathe + `_ (necessary if you want to build the + documentation (turned off by default) +- `Eigen `_ (3.3.0 or higher). If you do not + install this, it will be downloaded automatically at configuration time, so + this is not strictly necessary. The download can be slow, though, so we + recommend installing it on your system. - The CMake curses graphical user interface (``ccmake``). -*µ*\Spectre requires a relatively modern compiler as it makes heavy use of C++14 features. It has successfully been compiled and tested using the following compilers under Linux +*µ*\Spectre requires a relatively modern compiler as it makes heavy use of C++14 + features. It has successfully been compiled and tested using the following + compilers under Linux - gcc-7.2 - gcc-6.4 - gcc-5.4 - clang-6.0 - clang-5.0 - clang-4.0 and using clang-4.0 under MacOS. -It does *not* compile on Intel's most recent compiler, as it is still lacking some C++14 support. Work-arounds are planned, but will have to wait for someone to pick up the `task `_. +It does *not* compile on Intel's most recent compiler, as it is still lacking +some C++14 support. Work-arounds are planned, but will have to wait for someone +to pick up the `task `_. -To compile, create a build folder and configure the CMake project. If you do this in the folder you cloned in the previous step, it can look for instance like this: +To compile, create a build folder and configure the CMake project. If you do +this in the folder you cloned in the previous step, it can look for instance +like this: .. code-block:: sh $ mkdir build-release $ cd build-release $ ccmake .. -Then, set the build type to ``Release`` to produce optimised code. *µ*\Spectre makes heavy use of expression templates, so optimisation is paramount. (As an example, the performance difference between code compiled in ``Debug`` and ``Release`` is about a factor 40 in simple linear elasticity.) +Then, set the build type to ``Release`` to produce optimised code. *µ*\Spectre +makes heavy use of expression templates, so optimisation is paramount. (As an +example, the performance difference between code compiled in ``Debug`` and +``Release`` is about a factor 40 in simple linear elasticity.) Finally, compile the library and the tests by running .. code-block:: sh $ make -j .. warning:: - When using the ``-j`` option to compile, be aware that compiling *µ*\Spectre uses quite a bit of RAM. If your machine start swapping at compile time, reduce the number of parallel compilations + When using the ``-j`` option to compile, be aware that compiling *µ*\Spectre + uses quite a bit of RAM. If your machine start swapping at compile time, + reduce the number of parallel compilations Running *µ*\Spectre ******************* -The easiest and intended way of using *µ*\Spectre is through its Python bindings. The following simple example computes the response of a two-dimensional stretched periodic RVE cell. The cell consist of a soft matrix with a circular hard inclusion. +The easiest and intended way of using *µ*\Spectre is through its Python +bindings. The following simple example computes the response of a +two-dimensional stretched periodic RVE cell. The cell consist of a soft matrix +with a circular hard inclusion. .. literalinclude:: ../../../bin/tutorial_example.py :language: python -More examples both both python and c++ executables can be found in the ``/bin`` folder. +More examples both python and c++ executables can be found in the ``/bin`` +folder. Getting help ************ -*µ*\Spectre is in a very early stage of development and the documentation is currently spotty. Also, there is no FAQ page yet. If you run into trouble, please contact us on the `µSpectre chat room `_ and someone will answer as soon as possible. You can also check the API :ref:`reference`. +*µ*\Spectre is in a very early stage of development and the documentation is + currently spotty. Also, there is no FAQ page yet. If you run into trouble, + please contact us on the `µSpectre chat room `_ and + someone will answer as soon as possible. You can also check the API + :ref:`reference`. Reporting Bugs ************** -If you think you found a bug, you are probably right. Please report it! The preferred way is for you to create a task on `µSpectre's workboard `_ and assign it to user ``junge``. Include steps to reproduce the bug if possible. Someone will answer as soon as possible. +If you think you found a bug, you are probably right. Please report it! The +preferred way is for you to create a task on `µSpectre's workboard +`_ and assign it to user +``junge``. Include steps to reproduce the bug if possible. Someone will answer +as soon as possible. Contribute ********** -We welcome contributions both for new features and bug fixes. New features must be documented and have unit tests. Please submit contributions for review as `Arcanist revisions `_. More detailed guidelines for submissions will follow soonᵀᴹ. +We welcome contributions both for new features and bug fixes. New features must +be documented and have unit tests. Please submit contributions for review as +`Arcanist revisions +`_. More +detailed guidelines for submissions will follow soonᵀᴹ. diff --git a/doc/dev-docs/source/MaterialLinearElasticGeneric.rst b/doc/dev-docs/source/MaterialLinearElasticGeneric.rst index d07d2b5..ef62739 100644 --- a/doc/dev-docs/source/MaterialLinearElasticGeneric.rst +++ b/doc/dev-docs/source/MaterialLinearElasticGeneric.rst @@ -1,40 +1,51 @@ Generic Linear Elastic Material ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -The generic linear elastic material is implemented in the classes :cpp:class:`muSpectre::MaterialLinearElasticGeneric1` and :cpp:class:`muSpectre::MaterialLinearElasticGeneric2`, and is defined solely by the elastic stiffness tensor :math:`\mathbb{C}`, which has to be specified in `Voigt notation `_. In the case of :cpp:class:`muSpectre::MaterialLinearElasticGeneric2`, additionally, a per-pixel eigenstrain :math:`\bar{\boldsymbol{\varepsilon}}` can be supplied. The constitutive relation between the Cauchy stress :math:`\boldsymbol{\sigma}` and the small strain tensor :math:`\boldsymbol{\varepsilon}` is given by +The generic linear elastic material is implemented in the classes +:cpp:class:`MaterialLinearElasticGeneric1` and +:cpp:class:`MaterialLinearElasticGeneric2`, and is defined solely by +the elastic stiffness tensor :math:`\mathbb{C}`, which has to be specified in +`Voigt notation `_. In the case of +:cpp:class:`MaterialLinearElasticGeneric2`, additionally, a per-pixel +eigenstrain :math:`\bar{\boldsymbol{\varepsilon}}` can be supplied. The +constitutive relation between the Cauchy stress :math:`\boldsymbol{\sigma}` and +the small strain tensor :math:`\boldsymbol{\varepsilon}` is given by .. math:: :nowrap: \begin{align} \boldsymbol{\sigma} &= \mathbb{C}:\boldsymbol{\varepsilon}\\ \sigma_{ij} &= C_{ijkl}\,\varepsilon_{kl}, \quad\mbox{for the simple version, and}\\ \boldsymbol{\sigma} &= \mathbb{C}:\left(\boldsymbol{\varepsilon}-\bar{\boldsymbol{\varepsilon}}\right)\\ \sigma_{ij} &= C_{ijkl}\,\left(\varepsilon_{kl} - \bar\varepsilon_{kl}\right) \quad\mbox{ for the version with eigenstrain} \end{align} -This implementation is convenient, as it covers all possible linear elastic behaviours, but it is by far not as efficient as :cpp:class:`muSpectre::MaterialLinearElastic1` for isotropic linear elasticity. +This implementation is convenient, as it covers all possible linear elastic +behaviours, but it is by far not as efficient as +:cpp:class:`MaterialLinearElastic1` for isotropic linear elasticity. This law can be used in both small strain and finite strain calculations. -The following snippet shows how to use this law in python to implement isotropic linear elasticity: +The following snippet shows how to use this law in python to implement isotropic +linear elasticity: .. code-block:: python C = np.array([[2 * mu + lam, lam, lam, 0, 0, 0], [ lam, 2 * mu + lam, lam, 0, 0, 0], [ lam, lam, 2 * mu + lam, 0, 0, 0], [ 0, 0, 0, mu, 0, 0], [ 0, 0, 0, 0, mu, 0], [ 0, 0, 0, 0, 0, mu]]) eigenstrain = np.array([[ 0, .01], [.01, 0]]) mat1 = muSpectre.material.MaterialLinearElasticGeneric1_3d.make( cell, "material", C) mat1.add_pixel(pixel) mat2 = muSpectre.material.MaterialLinearElasticGeneric2_3d.make( cell, "material", C) mat2.add_pixel(pixel, eigenstrain) diff --git a/doc/dev-docs/source/NewMaterial.rst b/doc/dev-docs/source/NewMaterial.rst index 600a281..fc75006 100644 --- a/doc/dev-docs/source/NewMaterial.rst +++ b/doc/dev-docs/source/NewMaterial.rst @@ -1,206 +1,250 @@ Writing a New Constitutive Law ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The abstraction for a constitutive law in *µ*\Spectre** is the ``Material``, and new such materials must inherit from the class :cpp:class:`muSpectre::MaterialBase`. Most often, however, it will be most convenient to inherit from the derived class :cpp:class:`muSpectre::MaterialMuSpectre`, as it implements a lot of the machinery that is commonly used by constitutive laws. This section describes how to implement a new constitutive law with internal variables (sometimes also called state variables). The example material implemented here is :cpp:class:`MaterialTutorial`, an objective linear elastic law with a distribution of eigenstrains as internal variables. The constitutive law is defined by the relationship between material (or second Piola-Kirchhoff) stress :math:`\mathbf{S}` and Green-Lagrange strain :math:`\mathbf{E}` .. math:: :nowrap: \begin{align} \mathbf{S} &= \mathbb C:\left(\mathbf{E}-\mathbf{e}\right), \\ S_{ij} &= C_{ijkl}\left(E_{kl}-e_{kl}\right), \\ \end{align} where :math:`\mathbb C` is the elastic stiffness tensor and :math:`\mathbf e` is the local eigenstrain. Note that the implementation sketched out here is the most convenient to quickly get started with using *µ*\Spectre**, but not the most efficient one. For a most efficient implementation, refer to the implementation of :cpp:class:`muSpectre::MaterialLinearElastic2`. The :cpp:class:`muSpectre::MaterialMuSpectre` class *************************************************** -The class :cpp:class:`muSpectre::MaterialMuSpectre` is defined in ``material_muSpectre_base.hh`` and takes three template parameters; - -#. ``class Material`` is a `CRTP `_ and names the material inheriting from it. The reason for this construction is that we want to avoid virtual method calls from :cpp:class:`muSpectre::MaterialMuSpectre` to its derived classes. Rather, we want :cpp:class:`muSpectre::MaterialMuSpectre` to be able to call methods of the inheriting class directly without runtime overhead. -#. ``Dim_t DimS`` defines the number of spatial dimensions of the problem, i.e., whether we are dealing with a two- or three-dimensional grid of pixels/voxels. -#. ``Dim_t DimM`` defines the number of dimensions of our material description. This value will typically be the same as ``DimS``, but in cases like generalised plane strain, we can for instance have a three three-dimensional material response in a two-dimensional pixel grid. +The class :cpp:class:`muSpectre::MaterialMuSpectre` is defined in +``material_muSpectre_base.hh`` and takes three template parameters; + +#. ``class Material`` is a `CRTP + `_ and + names the material inheriting from it. The reason for this construction is + that we want to avoid virtual method calls from + :cpp:class:`muSpectre::MaterialMuSpectre` to its derived classes. Rather, we + want :cpp:class:`muSpectre::MaterialMuSpectre` to be able to call methods of + the inheriting class directly without runtime overhead. +#. ``Dim_t DimS`` defines the number of spatial dimensions of the problem, i.e., + whether we are dealing with a two- or three-dimensional grid of + pixels/voxels. +#. ``Dim_t DimM`` defines the number of dimensions of our material + description. This value will typically be the same as ``DimS``, but in cases + like generalised plane strain, we can for instance have a three + three-dimensional material response in a two-dimensional pixel grid. The main job of :cpp:class:`muSpectre::MaterialMuSpectre` is to -#. loop over all pixels to which this material has been assigned, transform the global gradient :math:`\mathbf{F}` (or small strain tensor :math:`\boldsymbol\varepsilon`) into the new material's required strain measure (e.g., the Green-Lagrange strain tensor :math:`\mathbf{E}`), -#. for each pixel evaluate the constitutive law by calling its ``evaluate_stress`` (computes the stress response) or ``evaluate_stress_tangent`` (computes both stress and consistent tangent) method with the local strain and internal variables, and finally -#. transform the stress (and possibly tangent) response from the material's stress measure into first Piola-Kirchhoff stress :math:`\mathbf{P}` (or Cauchy stress :math:`\boldsymbol\sigma` in small strain). - -In order to use these facilities, the new material needs to inherit from :cpp:class:`muSpectre::MaterialMuSpectre` (where we calculation of the response) and specialise the type :cpp:class:`muSpectre::MaterialMuSpectre_traits` (where we tell :cpp:class:`muSpectre::MaterialMuSpectre` how to use the new material). These two steps are described here for our example material. +#. loop over all pixels to which this material has been assigned, transform the + global gradient :math:`\mathbf{F}` (or small strain tensor + :math:`\boldsymbol\varepsilon`) into the new material's required strain + measure (e.g., the Green-Lagrange strain tensor :math:`\mathbf{E}`), +#. for each pixel evaluate the constitutive law by calling its + ``evaluate_stress`` (computes the stress response) or + ``evaluate_stress_tangent`` (computes both stress and consistent tangent) + method with the local strain and internal variables, and finally +#. transform the stress (and possibly tangent) response from the material's + stress measure into first Piola-Kirchhoff stress :math:`\mathbf{P}` (or + Cauchy stress :math:`\boldsymbol\sigma` in small strain). + +In order to use these facilities, the new material needs to inherit from +:cpp:class:`muSpectre::MaterialMuSpectre` (where we calculation of the response) +and specialise the type :cpp:class:`muSpectre::MaterialMuSpectre_traits` (where +we tell :cpp:class:`muSpectre::MaterialMuSpectre` how to use the new +material). These two steps are described here for our example material. Specialising the :cpp:class:`muSpectre::MaterialMuSpectre_traits` structure -*************************************************************************** -This structure is templated by the new material (in this case :cpp:class:`MaterialTutorial`) and needs to specify - -#. the types used to communicate per-pixel strains, stresses and stiffness tensors to the material (i.e., whether you want to get maps to `Eigen` matrices or raw pointers, or ...). Here we will use the convenient :cpp:type:`muSpectre::MatrixFieldMap` for strains and stresses, and :cpp:type:`muSpectre::T4MatrixFieldMap` for the stiffness. Look through the classes deriving from :cpp:type:`muSpectre::FieldMap` for all available options. -#. the strain measure that is expected (e.g., gradient, Green-Lagrange strain, left Cauchy-Green strain, etc.). Here we will use Green-Lagrange strain. The options are defined by the enum :cpp:enum:`muSpectre::StrainMeasure`. -#. the stress measure that is computed by the law (e.g., Cauchy, first Piola-Kirchhoff, etc,). Here, it will be first Piola-Kirchhoff stress. The available options are defined by the enum :cpp:enum:`muSpectre::StressMeasure`. +*************************************************************************** This +structure is templated by the new material (in this case +:cpp:class:`MaterialTutorial`) and needs to specify + +#. the types used to communicate per-pixel strains, stresses and stiffness + tensors to the material (i.e., whether you want to get maps to `Eigen` + matrices or raw pointers, or ...). Here we will use the convenient + :cpp:type:`muSpectre::MatrixFieldMap` for strains and stresses, and + :cpp:type:`muSpectre::T4MatrixFieldMap` for the stiffness. Look through the + classes deriving from :cpp:type:`muSpectre::FieldMap` for all available + options. +#. the strain measure that is expected (e.g., gradient, Green-Lagrange strain, + left Cauchy-Green strain, etc.). Here we will use Green-Lagrange strain. The + options are defined by the enum :cpp:enum:`muSpectre::StrainMeasure`. +#. the stress measure that is computed by the law (e.g., Cauchy, first + Piola-Kirchhoff, etc,). Here, it will be first Piola-Kirchhoff stress. The + available options are defined by the enum + :cpp:enum:`muSpectre::StressMeasure`. Our traits look like this (assuming we are in the namespace ``muSpectre``:: template struct MaterialMuSpectre_traits> { //! global field collection using GFieldCollection_t = typename GlobalFieldCollection; //! expected map type for strain fields using StrainMap_t = MatrixFieldMap; //! expected map type for stress fields using StressMap_t = MatrixFieldMap; //! expected map type for tangent stiffness fields using TangentMap_t = T4MatrixFieldMap; //! declare what type of strain measure your law takes as input constexpr static auto strain_measure{StrainMeasure::GreenLagrange}; //! declare what type of stress measure your law yields as output constexpr static auto stress_measure{StressMeasure::PK2}; //! local field_collections used for internals using LFieldColl_t = LocalFieldCollection; //! local strain type using LStrainMap_t = MatrixFieldMap; //! elasticity with eigenstrain using InternalVariables = std::tuple; }; Implementing the new material ***************************** -The new law needs to implement the methods ``add_pixel``, ``get_internals``, ``evaluate_stress``, and ``evaluate_stress_tangent``. Below is a commented example header:: +The new law needs to implement the methods ``add_pixel``, ``get_internals``, +``evaluate_stress``, and ``evaluate_stress_tangent``. Below is a commented +example header:: template class MaterialTutorial: public MaterialMuSpectre, DimS, DimM> { public: //! traits of this material using traits = MaterialMuSpectre_traits; //! Type of container used for storing eigenstrain using InternalVariables = typename traits::InternalVariables; //! Construct by name, Young's modulus and Poisson's ratio MaterialTutorial(std::string name, Real young, Real poisson); /** * evaluates second Piola-Kirchhoff stress given the Green-Lagrange * strain (or Cauchy stress if called with a small strain tensor) */ template inline decltype(auto) evaluate_stress(s_t && E, eigen_s_t && E_eig); /** * evaluates both second Piola-Kirchhoff stress and stiffness given * the Green-Lagrange strain (or Cauchy stress and stiffness if * called with a small strain tensor) */ template inline decltype(auto) evaluate_stress_tangent(s_t && E, eigen_s_t && E_eig); /** * return the internals tuple (needed by `muSpectre::MaterialMuSpectre`) */ InternalVariables & get_internals() { return this->internal_variables;}; /** * overload add_pixel to write into eigenstrain */ void add_pixel(const Ccoord_t & pixel, const Eigen::Matrix & E_eig); protected: //! stiffness tensor T4Mat C; //! storage for eigenstrain using Field_t = TensorField, Real, secondOrder, DimM>; Field_t & eigen_field; //!< field of eigenstrains //! tuple for iterable eigen_field InternalVariables internal_variables; private: }; A possible implementation for the constructor would be:: template MaterialTutorial::MaterialTutorial(std::string name, Real young, Real poisson) :MaterialMuSpectre(name) { // Lamé parameters Real lambda{young*poisson/((1+poisson)*(1-2*poisson))}; Real mu{young/(2*(1+poisson))}; // Kronecker delta Eigen::Matrix del{Eigen::Matrix::Identity()}; // fill the stiffness tensor this->C.setZero(); for (Dim_t i = 0; i < DimM; ++i) { for (Dim_t j = 0; j < DimM; ++j) { for (Dim_t k = 0; k < DimM; ++k) { for (Dim_t l = 0; l < DimM; ++l) { get(this->C, i, j, k, l) += (lambda * del(i,j)*del(k,l) + mu * (del(i,k)*del(j,l) + del(i,l)*del(j,k))); } } } } } -as an exercise, you could check how :cpp:class:`muSpectre::MaterialLinearElastic1` uses *µ*\Spectre**'s materials toolbox (in namespace ``MatTB``) to compute :math:`\mathbb C` in a much more convenient fashion. The evaluation of the stress could be (here, we make use of the ``Matrices`` namespace that defines common tensor algebra operations):: +as an exercise, you could check how +:cpp:class:`muSpectre::MaterialLinearElastic1` uses *µ*\Spectre**'s materials +toolbox (in namespace ``MatTB``) to compute :math:`\mathbb C` in a much more +convenient fashion. The evaluation of the stress could be (here, we make use of +the ``Matrices`` namespace that defines common tensor algebra operations):: template template decltype(auto) MaterialTutorial:: evaluate_stress(s_t && E, eigen_s_t && E_eig) { return Matrices::tens_mult(this->C, E-E_eig); } The remaining two methods are straight-forward:: template template decltype(auto) MaterialTutorial:: evaluate_stress_tangent(s_t && E, eigen_s_t && E_eig) { return return std::make_tuple (evaluate_stress(E, E_eig), this->C); } template InternalVariables & MaterialTutorial::get_internals() { return this->internal_variables; } -Note that the methods ``evaluate_stress`` and ``evaluate_stress_tangent`` need to be in the header, as both their input parameter types and output type depend on the compile-time context. +Note that the methods ``evaluate_stress`` and ``evaluate_stress_tangent`` need +to be in the header, as both their input parameter types and output type depend +on the compile-time context. diff --git a/language_bindings/python/bind_py_common.cc b/language_bindings/python/bind_py_common.cc index b143661..4289614 100644 --- a/language_bindings/python/bind_py_common.cc +++ b/language_bindings/python/bind_py_common.cc @@ -1,164 +1,174 @@ /** * @file bind_py_common.cc * * @author Till Junge * * @date 08 Jan 2018 * * @brief Python bindings for the common part of µSpectre * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "common/common.hh" #include "common/ccoord_operations.hh" #include #include #include using muSpectre::Dim_t; using muSpectre::Real; using muSpectre::StressMeasure; using muSpectre::StrainMeasure; using muSpectre::Formulation; using pybind11::literals::operator""_a; namespace py = pybind11; template void add_get_cube_helper(py::module & mod) { std::stringstream name{}; name << "get_" << dim << "d_cube"; mod.def(name.str().c_str(), &muSpectre::CcoordOps::get_cube, "size"_a, "return a Ccoord with the value 'size' repeated in each dimension"); } template void add_get_hermitian_helper(py::module & mod) { mod.def("get_hermitian_sizes", &muSpectre::CcoordOps::get_hermitian_sizes, "full_sizes"_a, "return the hermitian sizes corresponding to the true sizes"); } template void add_get_ccoord_helper(py::module & mod) { using Ccoord = muSpectre::Ccoord_t; mod.def( "get_domain_ccoord", [](Ccoord resolutions, Dim_t index) { return muSpectre::CcoordOps::get_ccoord(resolutions, Ccoord{}, index); }, "resolutions"_a, "i"_a, "return the cell coordinate corresponding to the i'th cell in a grid of " "shape resolutions"); } void add_get_cube(py::module & mod) { add_get_cube_helper(mod); add_get_cube_helper(mod); add_get_cube_helper(mod); add_get_cube_helper(mod); add_get_hermitian_helper(mod); add_get_hermitian_helper(mod); add_get_ccoord_helper(mod); add_get_ccoord_helper(mod); } template void add_get_index_helper(py::module & mod) { using Ccoord = muSpectre::Ccoord_t; mod.def("get_domain_index", [](Ccoord sizes, Ccoord ccoord) { return muSpectre::CcoordOps::get_index(sizes, Ccoord{}, ccoord); }, "sizes"_a, "ccoord"_a, "return the linear index corresponding to grid point 'ccoord' in a " "grid of size 'sizes'"); } void add_get_index(py::module & mod) { add_get_index_helper(mod); add_get_index_helper(mod); } template void add_Pixels_helper(py::module & mod) { std::stringstream name{}; name << "Pixels" << dim << "d"; using Ccoord = muSpectre::Ccoord_t; py::class_> Pixels(mod, name.str().c_str()); Pixels.def(py::init()); } void add_Pixels(py::module & mod) { add_Pixels_helper(mod); add_Pixels_helper(mod); } void add_common(py::module & mod) { py::enum_(mod, "Formulation") .value("finite_strain", Formulation::finite_strain) .value("small_strain", Formulation::small_strain); py::enum_(mod, "StressMeasure") .value("Cauchy", StressMeasure::Cauchy) .value("PK1", StressMeasure::PK1) .value("PK2", StressMeasure::PK2) .value("Kirchhoff", StressMeasure::Kirchhoff) .value("Biot", StressMeasure::Biot) .value("Mandel", StressMeasure::Mandel) .value("no_stress_", StressMeasure::no_stress_); py::enum_(mod, "StrainMeasure") .value("Gradient", StrainMeasure::Gradient) .value("Infinitesimal", StrainMeasure::Infinitesimal) .value("GreenLagrange", StrainMeasure::GreenLagrange) .value("Biot", StrainMeasure::Biot) .value("Log", StrainMeasure::Log) .value("Almansi", StrainMeasure::Almansi) .value("RCauchyGreen", StrainMeasure::RCauchyGreen) .value("LCauchyGreen", StrainMeasure::LCauchyGreen) .value("no_strain_", StrainMeasure::no_strain_); py::enum_(mod, "FFT_PlanFlags") .value("estimate", muSpectre::FFT_PlanFlags::estimate) .value("measure", muSpectre::FFT_PlanFlags::measure) .value("patient", muSpectre::FFT_PlanFlags::patient); + py::enum_( + mod, "FiniteDiff", + "Distinguishes between different options of numerical differentiation;\n " + " 1) 'forward' finite differences: ∂f/∂x ≈ (f(x+Δx) - f(x))/Δx\n 2) " + "'backward' finite differences: ∂f/∂x ≈ (f(x) - f(x-Δx))/Δx\n 3) " + "'centred' finite differences: ∂f/∂x ≈ (f(x+Δx) - f(x-Δx))/2Δx") + .value("forward", muSpectre::FiniteDiff::forward) + .value("backward", muSpectre::FiniteDiff::backward) + .value("centred", muSpectre::FiniteDiff::centred); + mod.def("banner", &muSpectre::banner, "name"_a, "year"_a, "copyright_holder"_a); add_get_cube(mod); add_Pixels(mod); add_get_index(mod); } diff --git a/language_bindings/python/bind_py_material.cc b/language_bindings/python/bind_py_material.cc index 5270398..bfeef87 100644 --- a/language_bindings/python/bind_py_material.cc +++ b/language_bindings/python/bind_py_material.cc @@ -1,159 +1,237 @@ /** * @file bind_py_material.cc * * @author Till Junge * * @date 09 Jan 2018 * * @brief python bindings for µSpectre's materials * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "common/common.hh" #include "materials/material_base.hh" +#include "materials/material_evaluator.hh" + #include #include #include #include #include using muSpectre::Dim_t; +using muSpectre::Real; using pybind11::literals::operator""_a; namespace py = pybind11; /* ---------------------------------------------------------------------- */ template void add_material_linear_elastic_generic1_helper(py::module & mod); /* ---------------------------------------------------------------------- */ template void add_material_linear_elastic_generic2_helper(py::module & mod); /** * python binding for the optionally objective form of Hooke's law */ -template +template void add_material_linear_elastic1_helper(py::module & mod); -template +template void add_material_linear_elastic2_helper(py::module & mod); -template +template void add_material_linear_elastic3_helper(py::module & mod); -template +template void add_material_linear_elastic4_helper(py::module & mod); template void add_material_hyper_elasto_plastic1_helper(py::module & mod); /* ---------------------------------------------------------------------- */ template class PyMaterialBase : public muSpectre::MaterialBase { public: /* Inherit the constructors */ using Parent = muSpectre::MaterialBase; using Parent::Parent; /* Trampoline (need one for each virtual function) */ void save_history_variables() override { PYBIND11_OVERLOAD_PURE(void, // Return type Parent, // Parent class save_history_variables); // Name of function in C++ // (must match Python name) } /* Trampoline (need one for each virtual function) */ void initialise() override { PYBIND11_OVERLOAD_PURE( void, // Return type Parent, // Parent class initialise); // Name of function in C++ (must match Python name) } void compute_stresses(const typename Parent::StrainField_t & F, typename Parent::StressField_t & P, muSpectre::Formulation form) override { PYBIND11_OVERLOAD_PURE( void, // Return type Parent, // Parent class compute_stresses, // Name of function in C++ (must match Python name) F, P, form); } void compute_stresses_tangent(const typename Parent::StrainField_t & F, typename Parent::StressField_t & P, typename Parent::TangentField_t & K, muSpectre::Formulation form) override { PYBIND11_OVERLOAD_PURE( void, /* Return type */ Parent, /* Parent class */ compute_stresses, /* Name of function in C++ (must match Python name) */ F, P, K, form); } }; -template +template +void add_material_evaluator(py::module & mod) { + std::stringstream name_stream{}; + name_stream << "MaterialEvaluator_" << Dim << "d"; + std::string name{name_stream.str()}; + + using MatEval_t = muSpectre::MaterialEvaluator; + py::class_(mod, name.c_str()) + .def(py::init>>()) + .def("save_history_variables", &MatEval_t::save_history_variables, + "for materials with state variables") + .def("evaluate_stress", + [](MatEval_t & mateval, py::EigenDRef & grad, + muSpectre::Formulation form) { + if ((grad.cols() != Dim) or (grad.rows() != Dim)) { + std::stringstream err{}; + err << "need matrix of shape (" << Dim << "×" << Dim + << ") but got (" << grad.rows() << "×" + << grad.cols() << ")."; + throw std::runtime_error(err.str()); + } + return mateval.evaluate_stress(grad, form); + }, + "strain"_a, "formulation"_a, + "Evaluates stress for a given strain and formulation " + "(Piola-Kirchhoff 1 stress as a function of the placement gradient " + "P = P(F) for formulation=Formulation.finite_strain and Cauchy " + "stress as a function of the infinitesimal strain tensor σ = σ(ε) " + "for formulation=Formulation.small_strain)") + .def("evaluate_stress_tangent", + [](MatEval_t & mateval, py::EigenDRef & grad, + muSpectre::Formulation form) { + if ((grad.cols() != Dim) or (grad.rows() != Dim)) { + std::stringstream err{}; + err << "need matrix of shape (" << Dim << "×" << Dim + << ") but got (" << grad.rows() << "×" + << grad.cols() << ")."; + throw std::runtime_error(err.str()); + } + return mateval.evaluate_stress_tangent(grad, form); + }, + "strain"_a, "formulation"_a, + "Evaluates stress and tangent moduli for a given strain and " + "formulation (Piola-Kirchhoff 1 stress as a function of the " + "placement gradient P = P(F) for " + "formulation=Formulation.finite_strain and Cauchy stress as a " + "function of the infinitesimal strain tensor σ = σ(ε) for " + "formulation=Formulation.small_strain). The tangent moduli are K = " + "∂P/∂F for formulation=Formulation.finite_strain and C = ∂σ/∂ε for " + "formulation=Formulation.small_strain.") + .def("estimate_tangent", + [](MatEval_t & evaluator, py::EigenDRef & grad, + muSpectre::Formulation form, const Real step, + const muSpectre::FiniteDiff diff_type) { + if ((grad.cols() != Dim) or (grad.rows() != Dim)) { + std::stringstream err{}; + err << "need matrix of shape (" << Dim << "×" << Dim + << ") but got (" << grad.rows() << "×" + << grad.cols() << ")."; + throw std::runtime_error(err.str()); + } + return evaluator.estimate_tangent(grad, form, step, diff_type); + }, + "strain"_a, "formulation"_a, "Delta_x"_a, + "difference_type"_a = muSpectre::FiniteDiff::centred, + "Numerical estimate of the tangent modulus using finite " + "differences. The finite difference scheme as well as the finite " + "step size can be chosen. If there are no special circumstances, " + "the default scheme of centred finite differences yields the most " + "accurate results at an increased computational cost."); +} + +template void add_material_helper(py::module & mod) { std::stringstream name_stream{}; - name_stream << "MaterialBase_" << dim << "d"; + name_stream << "MaterialBase_" << Dim << "d"; std::string name{name_stream.str()}; - using Material = muSpectre::MaterialBase; - using MaterialTrampoline = PyMaterialBase; - using FC_t = muSpectre::LocalFieldCollection; - using FCBase_t = muSpectre::FieldCollectionBase; + using Material = muSpectre::MaterialBase; + using MaterialTrampoline = PyMaterialBase; + using FC_t = muSpectre::LocalFieldCollection; + using FCBase_t = muSpectre::FieldCollectionBase; - py::class_(mod, + py::class_>(mod, name.c_str()) .def(py::init()) .def("save_history_variables", &Material::save_history_variables) .def("list_fields", &Material::list_fields) .def("get_real_field", &Material::get_real_field, "field_name"_a, py::return_value_policy::reference_internal) .def("size", &Material::size) .def("add_pixel", - [](Material & mat, muSpectre::Ccoord_t pix) { + [](Material & mat, muSpectre::Ccoord_t pix) { mat.add_pixel(pix); }, "pixel"_a) .def_property_readonly("collection", [](Material & material) -> FCBase_t & { return material.get_collection(); }, "returns the field collection containing internal " "fields of this material"); - add_material_linear_elastic1_helper(mod); - add_material_linear_elastic2_helper(mod); - add_material_linear_elastic3_helper(mod); - add_material_linear_elastic4_helper(mod); - add_material_hyper_elasto_plastic1_helper(mod); - add_material_linear_elastic_generic1_helper(mod); - add_material_linear_elastic_generic2_helper(mod); + add_material_linear_elastic1_helper(mod); + add_material_linear_elastic2_helper(mod); + add_material_linear_elastic3_helper(mod); + add_material_linear_elastic4_helper(mod); + add_material_hyper_elasto_plastic1_helper(mod); + add_material_linear_elastic_generic1_helper(mod); + add_material_linear_elastic_generic2_helper(mod); + + add_material_evaluator(mod); } void add_material(py::module & mod) { auto material{mod.def_submodule("material")}; material.doc() = "bindings for constitutive laws"; add_material_helper(material); add_material_helper(material); } diff --git a/language_bindings/python/bind_py_material_hyper_elasto_plastic1.cc b/language_bindings/python/bind_py_material_hyper_elasto_plastic1.cc index ae201fa..6f022e0 100644 --- a/language_bindings/python/bind_py_material_hyper_elasto_plastic1.cc +++ b/language_bindings/python/bind_py_material_hyper_elasto_plastic1.cc @@ -1,72 +1,78 @@ /** * @file bind_py_material_hyper_elasto_plastic1.cc * * @author Till Junge * * @date 29 Oct 2018 * * @brief python binding for MaterialHyperElastoPlastic1 * * Copyright © 2018 Till Junge * * µSpectre 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, or (at * your option) any later version. * * µSpectre 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 µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "common/common.hh" #include "materials/stress_transformations_Kirchhoff.hh" #include "materials/material_hyper_elasto_plastic1.hh" #include "cell/cell_base.hh" #include #include #include #include #include using muSpectre::Dim_t; using muSpectre::Real; using pybind11::literals::operator""_a; namespace py = pybind11; /** * python binding for the optionally objective form of Hooke's law * with per-pixel elastic properties */ template void add_material_hyper_elasto_plastic1_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "MaterialHyperElastoPlastic1_" << Dim << "d"; const auto name{name_stream.str()}; using Mat_t = muSpectre::MaterialHyperElastoPlastic1; using Cell_t = muSpectre::CellBase; - py::class_>(mod, name.c_str()) + py::class_, std::shared_ptr>( + mod, name.c_str()) .def_static("make", [](Cell_t & cell, std::string name, Real Young, Real Poisson, Real tau_y0, Real h) -> Mat_t & { return Mat_t::make(cell, name, Young, Poisson, tau_y0, h); }, "cell"_a, "name"_a, "YoungModulus"_a, "PoissonRatio"_a, "τ_y₀"_a, "h"_a, py::return_value_policy::reference, - py::keep_alive<1, 0>()); + py::keep_alive<1, 0>()) + .def_static("make_evaluator", + [](Real Young, Real Poisson, Real tau_y0, Real h) { + return Mat_t::make_evaluator(Young, Poisson, tau_y0, h); + }, + "YoungModulus"_a, "PoissonRatio"_a, "τ_y₀"_a, "h"_a); } template void add_material_hyper_elasto_plastic1_helper(py::module &); template void add_material_hyper_elasto_plastic1_helper(py::module &); diff --git a/language_bindings/python/bind_py_material_linear_elastic1.cc b/language_bindings/python/bind_py_material_linear_elastic1.cc index 4f40993..7a6996b 100644 --- a/language_bindings/python/bind_py_material_linear_elastic1.cc +++ b/language_bindings/python/bind_py_material_linear_elastic1.cc @@ -1,68 +1,74 @@ /** * @file bind_py_material_linear_elastic1.cc * * @author Till Junge * * @date 29 Oct 2018 * * @brief python bindings for MaterialLinearElastic1 * * Copyright © 2018 Till Junge * * µSpectre 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, or (at * your option) any later version. * * µSpectre 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 µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "common/common.hh" #include "materials/material_linear_elastic1.hh" #include "cell/cell_base.hh" #include #include #include #include #include using muSpectre::Dim_t; using muSpectre::Real; using pybind11::literals::operator""_a; namespace py = pybind11; /** * python binding for the optionally objective form of Hooke's law */ template void add_material_linear_elastic1_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "MaterialLinearElastic1_" << dim << 'd'; const auto name{name_stream.str()}; using Mat_t = muSpectre::MaterialLinearElastic1; using Sys_t = muSpectre::CellBase; - py::class_>(mod, name.c_str()) + py::class_, std::shared_ptr>( + mod, name.c_str()) .def_static("make", [](Sys_t & sys, std::string n, Real e, Real p) -> Mat_t & { return Mat_t::make(sys, n, e, p); }, "cell"_a, "name"_a, "Young"_a, "Poisson"_a, - py::return_value_policy::reference, py::keep_alive<1, 0>()); + py::return_value_policy::reference, py::keep_alive<1, 0>()) + .def_static("make_evaluator", + [](Real e, Real p) { + return Mat_t::make_evaluator(e, p); + }, + "Young"_a, "Poisson"_a); } template void add_material_linear_elastic1_helper(py::module &); template void add_material_linear_elastic1_helper(py::module &); diff --git a/language_bindings/python/bind_py_material_linear_elastic2.cc b/language_bindings/python/bind_py_material_linear_elastic2.cc index 9ac60ce..a1fb000 100644 --- a/language_bindings/python/bind_py_material_linear_elastic2.cc +++ b/language_bindings/python/bind_py_material_linear_elastic2.cc @@ -1,77 +1,83 @@ /** * @file bind_py_material_linear_elastic2.cc * * @author Till Junge * * @date 29 Oct 2018 * * @brief python bindings for MaterialLinearElastic2 * * Copyright © 2018 Till Junge * * µSpectre 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, or (at * your option) any later version. * * µSpectre 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 µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "common/common.hh" #include "materials/material_linear_elastic2.hh" #include "cell/cell_base.hh" #include #include #include #include #include using muSpectre::Dim_t; using muSpectre::Real; namespace py = pybind11; using pybind11::literals::operator""_a; /** * python binding for the optionally objective form of Hooke's law * with a per pixel eigenstrain */ template void add_material_linear_elastic2_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "MaterialLinearElastic2_" << dim << 'd'; const auto name{name_stream.str()}; using Mat_t = muSpectre::MaterialLinearElastic2; using Sys_t = muSpectre::CellBase; - py::class_>(mod, name.c_str()) + py::class_, std::shared_ptr>( + mod, name.c_str()) .def_static("make", [](Sys_t & sys, std::string n, Real e, Real p) -> Mat_t & { return Mat_t::make(sys, n, e, p); }, "cell"_a, "name"_a, "Young"_a, "Poisson"_a, py::return_value_policy::reference, py::keep_alive<1, 0>()) .def("add_pixel", [](Mat_t & mat, muSpectre::Ccoord_t pix, py::EigenDRef & eig) { Eigen::Matrix eig_strain{eig}; mat.add_pixel(pix, eig_strain); }, - "pixel"_a, "eigenstrain"_a); + "pixel"_a, "eigenstrain"_a) + .def_static("make_evaluator", + [](Real e, Real p) { + return Mat_t::make_evaluator(e, p); + }, + "Young"_a, "Poisson"_a); } template void add_material_linear_elastic2_helper(py::module &); template void add_material_linear_elastic2_helper(py::module &); diff --git a/language_bindings/python/bind_py_material_linear_elastic3.cc b/language_bindings/python/bind_py_material_linear_elastic3.cc index 8778196..d11346b 100644 --- a/language_bindings/python/bind_py_material_linear_elastic3.cc +++ b/language_bindings/python/bind_py_material_linear_elastic3.cc @@ -1,74 +1,77 @@ /** * @file bind_py_material_linear_elastic3.cc * * @author Till Junge * * @date 29 Oct 2018 * * @brief python bindings for MaterialLinearElastic3 * * Copyright © 2018 Till Junge * * µSpectre 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, or (at * your option) any later version. * * µSpectre 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 µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "common/common.hh" #include "materials/material_linear_elastic3.hh" #include "cell/cell_base.hh" #include #include #include #include #include using muSpectre::Dim_t; using muSpectre::Real; using pybind11::literals::operator""_a; namespace py = pybind11; /** * python binding for the optionally objective form of Hooke's law * with per-pixel elastic properties */ template void add_material_linear_elastic3_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "MaterialLinearElastic3_" << dim << 'd'; const auto name{name_stream.str()}; using Mat_t = muSpectre::MaterialLinearElastic3; using Sys_t = muSpectre::CellBase; - py::class_>(mod, name.c_str()) + py::class_, std::shared_ptr>( + mod, name.c_str()) .def(py::init(), "name"_a) .def_static("make", [](Sys_t & sys, std::string n) -> Mat_t & { return Mat_t::make(sys, n); }, "cell"_a, "name"_a, py::return_value_policy::reference, py::keep_alive<1, 0>()) .def("add_pixel", [](Mat_t & mat, muSpectre::Ccoord_t pix, Real Young, Real Poisson) { mat.add_pixel(pix, Young, Poisson); }, - "pixel"_a, "Young"_a, "Poisson"_a); + "pixel"_a, "Young"_a, "Poisson"_a) + .def_static("make_evaluator", + []() { return Mat_t::make_evaluator(); }); } template void add_material_linear_elastic3_helper(py::module &); template void add_material_linear_elastic3_helper(py::module &); diff --git a/language_bindings/python/bind_py_material_linear_elastic4.cc b/language_bindings/python/bind_py_material_linear_elastic4.cc index fa3f782..6134460 100644 --- a/language_bindings/python/bind_py_material_linear_elastic4.cc +++ b/language_bindings/python/bind_py_material_linear_elastic4.cc @@ -1,74 +1,77 @@ /** * @file bind_py_material_linear_elastic4.cc * * @author Till Junge * * @date 29 Oct 2018 * * @brief python binding for MaterialLinearElastic4 * * Copyright © 2018 Till Junge * * µSpectre 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, or (at * your option) any later version. * * µSpectre 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 µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. */ #include "common/common.hh" #include "materials/material_linear_elastic4.hh" #include "cell/cell_base.hh" #include #include #include #include #include using muSpectre::Dim_t; using muSpectre::Real; using pybind11::literals::operator""_a; namespace py = pybind11; /** * python binding for the optionally objective form of Hooke's law * with per-pixel elastic properties */ template void add_material_linear_elastic4_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "MaterialLinearElastic4_" << dim << 'd'; const auto name{name_stream.str()}; using Mat_t = muSpectre::MaterialLinearElastic4; using Sys_t = muSpectre::CellBase; - py::class_>(mod, name.c_str()) + py::class_, std::shared_ptr>( + mod, name.c_str()) .def(py::init(), "name"_a) .def_static("make", [](Sys_t & sys, std::string n) -> Mat_t & { return Mat_t::make(sys, n); }, "cell"_a, "name"_a, py::return_value_policy::reference, py::keep_alive<1, 0>()) .def("add_pixel", [](Mat_t & mat, muSpectre::Ccoord_t pix, Real Young, Real Poisson) { mat.add_pixel(pix, Young, Poisson); }, - "pixel"_a, "Young"_a, "Poisson"_a); + "pixel"_a, "Young"_a, "Poisson"_a) + .def_static("make_evaluator", + []() { return Mat_t::make_evaluator(); }); } template void add_material_linear_elastic4_helper(py::module &); template void add_material_linear_elastic4_helper(py::module &); diff --git a/language_bindings/python/bind_py_material_linear_elastic_generic.cc b/language_bindings/python/bind_py_material_linear_elastic_generic.cc index 5a55bde..9ca6771 100644 --- a/language_bindings/python/bind_py_material_linear_elastic_generic.cc +++ b/language_bindings/python/bind_py_material_linear_elastic_generic.cc @@ -1,134 +1,150 @@ /** * @file bind_py_material_linear_elastic_generic.cc * * @author Till Junge * * @date 20 Dec 2018 * * @brief bindings for the generic linear elastic law defined by its stiffness * tensor * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "materials/material_linear_elastic_generic1.hh" #include "materials/material_linear_elastic_generic2.hh" #include "cell/cell_base.hh" #include #include #include #include #include using namespace muSpectre; // NOLINT // TODO(junge): figure this out namespace py = pybind11; using namespace pybind11::literals; // NOLINT: recommended usage /** * python binding for the generic linear elastic material */ template void add_material_linear_elastic_generic1_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "MaterialLinearElasticGeneric1_" << Dim << "d"; const auto name{name_stream.str()}; using Mat_t = MaterialLinearElasticGeneric1; using Cell_t = CellBase; - py::class_>(mod, name.c_str()) + py::class_, std::shared_ptr>( + mod, name.c_str()) .def_static( "make", [](Cell_t & cell, std::string name, const py::EigenDRef< Eigen::Matrix> & elastic_tensor) -> Mat_t & { return Mat_t::make(cell, name, elastic_tensor); }, "cell"_a, "name"_a, "elastic_tensor"_a, py::return_value_policy::reference, py::keep_alive<1, 0>(), "Factory function returning a MaterialLinearElastic instance. " "The elastic tensor has to be specified in Voigt notation.") .def("add_pixel", [](Mat_t & mat, Ccoord_t pix) { mat.add_pixel(pix); }, "pixel"_a, "Register a new pixel to this material. subsequent evaluations of " "the " "stress and tangent in the cell will use this constitutive law for " "this " "particular pixel") - .def("size", &Mat_t::size); + .def("size", &Mat_t::size) + .def_static("make_evaluator", + [](const py::EigenDRef< + Eigen::Matrix> & + elastic_tensor) { + return Mat_t::make_evaluator(elastic_tensor); + }, + "elastic_tensor"_a); } /** * python binding for the generic linear elastic material with eigenstcain */ template void add_material_linear_elastic_generic2_helper(py::module & mod) { std::stringstream name_stream{}; name_stream << "MaterialLinearElasticGeneric2_" << Dim << "d"; const auto name{name_stream.str()}; using Mat_t = MaterialLinearElasticGeneric2; using Cell_t = CellBase; - py::class_>(mod, name.c_str()) + py::class_, std::shared_ptr>( + mod, name.c_str()) .def_static( "make", [](Cell_t & cell, std::string name, const py::EigenDRef< Eigen::Matrix> & elastic_tensor) -> Mat_t & { return Mat_t::make(cell, name, elastic_tensor); }, "cell"_a, "name"_a, "elastic_tensor"_a, py::return_value_policy::reference, py::keep_alive<1, 0>(), "Factory function returning a MaterialLinearElastic instance. " "The elastic tensor has to be specified in Voigt notation.") .def("add_pixel", [](Mat_t & mat, Ccoord_t pix) { mat.add_pixel(pix); }, "pixel"_a, "Register a new pixel to this material. Subsequent evaluations of " "the stress and tangent in the cell will use this constitutive law " "for this particular pixel") .def("add_pixel", [](Mat_t & mat, Ccoord_t pix, py::EigenDRef & eig) { Eigen::Matrix eig_strain{eig}; mat.add_pixel(pix, eig_strain); }, "pixel"_a, "eigenstrain"_a, "Register a new pixel to this material and assign the eigenstrain. " "Subsequent Evaluations of the stress and tangent in the cell will " "use this constitutive law for this particular pixel") - .def("size", &Mat_t::size); + .def("size", &Mat_t::size) + .def_static("make_evaluator", + [](const py::EigenDRef< + Eigen::Matrix> & + elastic_tensor) { + return Mat_t::make_evaluator(elastic_tensor); + }, + "elastic_tensor"_a); } template void add_material_linear_elastic_generic1_helper(py::module &); template void add_material_linear_elastic_generic1_helper(py::module &); template void add_material_linear_elastic_generic2_helper(py::module &); template void add_material_linear_elastic_generic2_helper(py::module &); diff --git a/language_bindings/python/muSpectre/__init__.py b/language_bindings/python/muSpectre/__init__.py index 7c0da4d..e6d39bc 100644 --- a/language_bindings/python/muSpectre/__init__.py +++ b/language_bindings/python/muSpectre/__init__.py @@ -1,104 +1,104 @@ #!/usr/bin/env python3 # -*- coding:utf-8 -*- """ @file __init__.py @author Lars Pastewka @date 21 Mar 2018 @brief Main entry point for muSpectre Python module Copyright © 2018 Till Junge µSpectre is free software; you can redistribute it and/or modify it under the terms of the GNU General Lesser Public License as published by the Free Software Foundation, either version 3, or (at your option) any later version. µSpectre 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 Lesser General Public License along with µSpectre; see the file COPYING. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. Additional permission under GNU GPL version 3 section 7 If you modify this Program, or any covered work, by linking or combining it with proprietary FFT implementations or numerical libraries, containing parts covered by the terms of those libraries' licenses, the licensors of this Program grant you additional permission to convey the resulting work. """ try: from mpi4py import MPI except ImportError: MPI = None import _muSpectre from _muSpectre import (Formulation, get_domain_ccoord, get_domain_index, get_hermitian_sizes, material, solvers, - FFT_PlanFlags) + FFT_PlanFlags, FiniteDiff) import muSpectre.fft _factories = {'fftw': ('CellFactory', False), 'fftwmpi': ('FFTWMPICellFactory', True), 'pfft': ('PFFTCellFactory', True), 'p3dfft': ('P3DFFTCellFactory', True)} def Cell(resolutions, lengths, formulation=Formulation.finite_strain, fft='fftw', communicator=None): """ Instantiate a muSpectre Cell class. Parameters ---------- resolutions: list Grid resolutions in the Cartesian directions. lengths: list Physical size of the cell in the Cartesian directions. formulation: Formulation Formulation for strains and stresses used by the solver. Options are `Formulation.finite_strain` and `Formulation.small_strain`. Finite strain formulation is the default. fft: string FFT engine to use. Options are 'fftw', 'fftwmpi', 'pfft' and 'p3dfft'. Default is 'fftw'. communicator: mpi4py communicator mpi4py communicator object passed to parallel FFT engines. Note that the default 'fftw' engine does not support parallel execution. Returns ------- cell: object Return a muSpectre Cell object. """ try: factory_name, is_parallel = _factories[fft] except KeyError: raise KeyError("Unknown FFT engine '{}'.".format(fft)) try: factory = _muSpectre.__dict__[factory_name] except KeyError: raise KeyError("FFT engine '{}' has not been compiled into the " "muSpectre library.".format(fft)) if is_parallel: if MPI is None: raise RuntimeError('Parallel solver requested but mpi4py could' ' not be imported.') if communicator is None: communicator = MPI.COMM_SELF return factory(resolutions, lengths, formulation, MPI._handleof(communicator)) else: if communicator is not None: raise ValueError("FFT engine '{}' does not support parallel " "execution.".format(fft)) return factory(resolutions, lengths, formulation) diff --git a/src/cell/cell_base.hh b/src/cell/cell_base.hh index f72f350..465d26e 100644 --- a/src/cell/cell_base.hh +++ b/src/cell/cell_base.hh @@ -1,679 +1,676 @@ /** * @file cell_base.hh * * @author Till Junge * * @date 01 Nov 2017 * * @brief Base class representing a unit cell cell with single * projection operator * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #ifndef SRC_CELL_CELL_BASE_HH_ #define SRC_CELL_CELL_BASE_HH_ #include "common/common.hh" #include "common/ccoord_operations.hh" #include "common/field.hh" #include "common/utilities.hh" #include "materials/material_base.hh" #include "fft/projection_base.hh" #include "cell/cell_traits.hh" #include #include #include #include #include namespace muSpectre { /** * Cell adaptors implement the matrix-vector multiplication and * allow the system to be used like a sparse matrix in * conjugate-gradient-type solvers */ - template - class CellAdaptor; + template class CellAdaptor; /** * Base class for cells that is not templated and therefore can be * in solvers that see cells as runtime-polymorphic objects. This * allows the use of standard * (i.e. spectral-method-implementation-agnostic) solvers, as for * instance the scipy solvers */ class Cell { public: //! sparse matrix emulation using Adaptor = CellAdaptor; //! dynamic vector type for interactions with numpy/scipy/solvers etc. using Vector_t = Eigen::Matrix; //! dynamic matrix type for setting strains using Matrix_t = Eigen::Matrix; //! dynamic generic array type for interaction with numpy, i/o, etc template using Array_t = Eigen::Array; //! ref to dynamic generic array - template - using Array_ref = Eigen::Map>; + template using Array_ref = Eigen::Map>; //! ref to constant vector using ConstVector_ref = Eigen::Map; //! output vector reference for solvers using Vector_ref = Eigen::Map; //! Default constructor Cell() = default; //! Copy constructor Cell(const Cell & other) = default; //! Move constructor Cell(Cell && other) = default; //! Destructor virtual ~Cell() = default; //! Copy assignment operator Cell & operator=(const Cell & other) = default; //! Move assignment operator Cell & operator=(Cell && other) = default; //! for handling double initialisations right bool is_initialised() const { return this->initialised; } //! returns the number of degrees of freedom in the cell virtual Dim_t get_nb_dof() const = 0; //! number of pixels in the cell virtual size_t size() const = 0; //! return the communicator object virtual const Communicator & get_communicator() const = 0; /** * formulation is hard set by the choice of the projection class */ virtual const Formulation & get_formulation() const = 0; /** * returns the material dimension of the problem */ virtual Dim_t get_material_dim() const = 0; /** * returns the number of rows and cols for the strain matrix type * (for full storage, the strain is stored in material_dim × * material_dim matrices, but in symmetriy storage, it is a column * vector) */ virtual std::array get_strain_shape() const = 0; /** * returns a writable map onto the strain field of this cell. This * corresponds to the unknowns in a typical solve cycle. */ virtual Vector_ref get_strain_vector() = 0; /** * returns a read-only map onto the stress field of this * cell. This corresponds to the intermediate (and finally, total) * solution in a typical solve cycle */ virtual ConstVector_ref get_stress_vector() const = 0; /** * evaluates and returns the stress for the currently set strain */ virtual ConstVector_ref evaluate_stress() = 0; /** * evaluates and returns the stress and stiffness for the currently set * strain */ virtual std::array evaluate_stress_tangent() = 0; /** * applies the projection operator in-place on the input vector */ virtual void apply_projection(Eigen::Ref vec) = 0; /** * evaluates the projection of the input field (this corresponds * to G:P in de Geus 2017, * http://dx.doi.org/10.1016/j.cma.2016.12.032). The first time, * this allocates the memory for the return value, and reuses it * on subsequent calls */ virtual Vector_ref evaluate_projection(Eigen::Ref P) = 0; /** * freezes all the history variables of the materials */ virtual void save_history_variables() = 0; /** * evaluates the directional and projected stiffness (this * corresponds to G:K:δF in de Geus 2017, * http://dx.doi.org/10.1016/j.cma.2016.12.032). It seems that * this operation needs to be implemented with a copy in oder to * be compatible with scipy and EigenCG etc (At the very least, * the copy is only made once) */ virtual Vector_ref evaluate_projected_directional_stiffness( Eigen::Ref delF) = 0; /** * returns a ref to a field named 'unique_name" of real values * managed by the cell. If the field does not yet exist, it is * created. * * @param unique_name name of the field. If the field already * exists, an array ref mapped onto it is returned. Else, a new * field with that name is created and returned- * * @param nb_components number of components to be stored *per * pixel*. For new fields any positive number can be chosen. When * accessing an existing field, this must correspond to the * existing field size, and a `std::runtime_error` is thrown if * this is not satisfied */ virtual Array_ref get_managed_real_array(std::string unique_name, size_t nb_components) = 0; /** * Convenience function to copy local (internal) fields of * materials into a global field. At least one of the materials in * the cell needs to contain an internal field named * `unique_name`. If multiple materials contain such a field, they * all need to be of same scalar type and same number of * components. This does not work for split pixel cells or * laminate pixel cells, as they can have multiple entries for the * same pixel. Pixels for which no field named `unique_name` * exists get an array of zeros. * * @param unique_name fieldname to fill the global field with. At * least one material must have such a field, or a * `std::runtime_error` is thrown */ virtual Array_ref get_globalised_internal_real_array(const std::string & unique_name) = 0; /** * Convenience function to copy local (internal) state fields * (current state) of materials into a global field. At least one * of the materials in the cell needs to contain an internal field * named `unique_name`. If multiple materials contain such a * field, they all need to be of same scalar type and same number * of components. This does not work for split pixel cells or * laminate pixel cells, as they can have multiple entries for the * same pixel. Pixels for which no field named `unique_name` * exists get an array of zeros. * * @param unique_name fieldname to fill the global field with. At * least one material must have such a field, or a * `std::runtime_error` is thrown */ virtual Array_ref get_globalised_current_real_array(const std::string & unique_name) = 0; /** * Convenience function to copy local (internal) state fields * (old state) of materials into a global field. At least one * of the materials in the cell needs to contain an internal field * named `unique_name`. If multiple materials contain such a * field, they all need to be of same scalar type and same number * of components. This does not work for split pixel cells or * laminate pixel cells, as they can have multiple entries for the * same pixel. Pixels for which no field named `unique_name` * exists get an array of zeros. * - * @param unique_name fieldname to fill the global field with. At - * least one material must have such a field, or a - * `std::runtime_error` is thrown + * @param unique_name fieldname to fill the global field with. At least one + * material must have such a field, or a `std::runtime_error` is thrown + * @param nb_steps_ago for history variables which remember more than a + * single previous value, `nb_steps_ago` can be used to specify which old + * value to access. */ virtual Array_ref get_globalised_old_real_array(const std::string & unique_name, int nb_steps_ago = 1) = 0; /** * set uniform strain (typically used to initialise problems */ virtual void set_uniform_strain(const Eigen::Ref &) = 0; //! get a sparse matrix view on the cell virtual Adaptor get_adaptor() = 0; protected: bool initialised{false}; //!< to handle double initialisation right private: }; //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) template class CellBase : public Cell { public: using Parent = Cell; using Ccoord = Ccoord_t; //!< cell coordinates type using Rcoord = Rcoord_t; //!< physical coordinates type //! global field collection using FieldCollection_t = GlobalFieldCollection; //! the collection is handled in a `std::unique_ptr` using Collection_ptr = std::unique_ptr; //! polymorphic base material type using Material_t = MaterialBase; //! materials handled through `std::unique_ptr`s using Material_ptr = std::unique_ptr; //! polymorphic base projection type using Projection_t = ProjectionBase; //! projections handled through `std::unique_ptr`s using Projection_ptr = std::unique_ptr; //! dynamic global fields - template - using Field_t = TypedField; + template using Field_t = TypedField; //! expected type for strain fields using StrainField_t = TensorField; //! expected type for stress fields using StressField_t = TensorField; //! expected type for tangent stiffness fields using TangentField_t = TensorField; //! combined stress and tangent field using FullResponse_t = std::tuple; //! iterator type over all cell pixel's using iterator = typename CcoordOps::Pixels::iterator; //! dynamic vector type for interactions with numpy/scipy/solvers etc. using Vector_t = typename Parent::Vector_t; //! ref to constant vector using ConstVector_ref = typename Parent::ConstVector_ref; //! output vector reference for solvers using Vector_ref = typename Parent::Vector_ref; //! dynamic array type for interactions with numpy/scipy/solvers, etc. - template - using Array_t = typename Parent::Array_t; + template using Array_t = typename Parent::Array_t; //! dynamic array type for interactions with numpy/scipy/solvers, etc. - template - using Array_ref = typename Parent::Array_ref; + template using Array_ref = typename Parent::Array_ref; //! sparse matrix emulation using Adaptor = Parent::Adaptor; //! Default constructor CellBase() = delete; //! constructor using sizes and resolution explicit CellBase(Projection_ptr projection); //! Copy constructor CellBase(const CellBase & other) = delete; //! Move constructor CellBase(CellBase && other); //! Destructor virtual ~CellBase() = default; //! Copy assignment operator CellBase & operator=(const CellBase & other) = delete; //! Move assignment operator CellBase & operator=(CellBase && other) = default; /** * Materials can only be moved. This is to assure exclusive * ownership of any material by this cell */ Material_t & add_material(Material_ptr mat); /** * returns a writable map onto the strain field of this cell. This * corresponds to the unknowns in a typical solve cycle. */ Vector_ref get_strain_vector() override; /** * returns a read-only map onto the stress field of this * cell. This corresponds to the intermediate (and finally, total) * solution in a typical solve cycle */ ConstVector_ref get_stress_vector() const override; /** * evaluates and returns the stress for the currently set strain */ ConstVector_ref evaluate_stress() override; /** * evaluates and returns the stress and stiffness for the currently set * strain */ std::array evaluate_stress_tangent() override; /** * evaluates the projection of the input field (this corresponds * do G:P in de Geus 2017, * http://dx.doi.org/10.1016/j.cma.2016.12.032). The first time, * this allocates the memory for the return value, and reuses it * on subsequent calls */ Vector_ref evaluate_projection(Eigen::Ref P) override; /** * evaluates the directional and projected stiffness (this * corresponds to G:K:δF in de Geus 2017, * http://dx.doi.org/10.1016/j.cma.2016.12.032). It seems that * this operation needs to be implemented with a copy in oder to * be compatible with scipy and EigenCG etc. (At the very least, * the copy is only made once) */ Vector_ref evaluate_projected_directional_stiffness( Eigen::Ref delF) override; //! return the template param DimM (required for polymorphic use of `Cell` Dim_t get_material_dim() const final { return DimM; } /** * returns the number of rows and cols for the strain matrix type * (for full storage, the strain is stored in material_dim × * material_dim matrices, but in symmetriy storage, it is a column * vector) */ std::array get_strain_shape() const final; /** * applies the projection operator in-place on the input vector */ void apply_projection(Eigen::Ref vec) final; /** * set uniform strain (typically used to initialise problems */ void set_uniform_strain(const Eigen::Ref &) override; /** * evaluates all materials */ FullResponse_t evaluate_stress_tangent(StrainField_t & F); /** * evaluate directional stiffness (i.e. G:K:δF or G:K:δε) */ StressField_t & directional_stiffness(const TangentField_t & K, const StrainField_t & delF, StressField_t & delP); /** * vectorized version for eigen solvers, no copy, but only works * when fields have ArrayStore=false */ Vector_ref directional_stiffness_vec(const Eigen::Ref & delF); /** * Evaluate directional stiffness into a temporary array and * return a copy. This is a costly and wasteful interface to * directional_stiffness and should only be used for debugging or * in the python interface */ Eigen::ArrayXXd directional_stiffness_with_copy(Eigen::Ref delF); /** * Convenience function circumventing the neeed to use the * underlying projection */ StressField_t & project(StressField_t & field); //! returns a ref to the cell's strain field StrainField_t & get_strain(); //! returns a ref to the cell's stress field const StressField_t & get_stress() const; //! returns a ref to the cell's tangent stiffness field const TangentField_t & get_tangent(bool create = false); //! returns a ref to a temporary field managed by the cell StrainField_t & get_managed_T2_field(std::string unique_name); //! returns a ref to a temporary field of real values managed by the cell Field_t & get_managed_real_field(std::string unique_name, size_t nb_components); /** * returns a Array ref to a temporary field of real values managed by the * cell */ Array_ref get_managed_real_array(std::string unique_name, size_t nb_components) final; /** * returns a global field filled from local (internal) fields of * the materials. see `Cell::get_globalised_internal_real_array` for * details. */ Field_t & get_globalised_internal_real_field(const std::string & unique_name); /** * returns a global field filled from local (internal) statefields of * the materials. see `Cell::get_globalised_current_real_array` for * details. */ Field_t & get_globalised_current_real_field(const std::string & unique_name); /** * returns a global field filled from local (internal) statefields of * the materials. see `Cell::get_globalised_old_real_array` for * details. */ Field_t & get_globalised_old_real_field(const std::string & unique_name, int nb_steps_ago = 1); //! see `Cell::get_globalised_internal_real_array` for details Array_ref get_globalised_internal_real_array( const std::string & unique_name) final; //! see `Cell::get_globalised_current_reald_array` for details Array_ref get_globalised_current_real_array( const std::string & unique_name) final; //! see `Cell::get_globalised_old_real_array` for details Array_ref get_globalised_old_real_array(const std::string & unique_name, int nb_steps_ago = 1) final; /** * general initialisation; initialises the projection and * fft_engine (i.e. infrastructure) but not the materials. These * need to be initialised separately */ void initialise(FFT_PlanFlags flags = FFT_PlanFlags::estimate); /** * for materials with state variables, these typically need to be * saved/updated an the end of each load increment, this function * calls this update for each material in the cell */ void save_history_variables() final; iterator begin(); //!< iterator to the first pixel iterator end(); //!< iterator past the last pixel //! number of pixels in the cell size_t size() const final { return pixels.size(); } //! return the subdomain resolutions of the cell const Ccoord & get_subdomain_resolutions() const { return this->subdomain_resolutions; } //! return the subdomain locations of the cell const Ccoord & get_subdomain_locations() const { return this->subdomain_locations; } //! return the domain resolutions of the cell const Ccoord & get_domain_resolutions() const { return this->domain_resolutions; } //! return the sizes of the cell const Rcoord & get_domain_lengths() const { return this->domain_lengths; } /** * formulation is hard set by the choice of the projection class */ const Formulation & get_formulation() const final { return this->projection->get_formulation(); } /** * get a reference to the projection object. should only be * required for debugging */ Eigen::Map get_projection() { return this->projection->get_operator(); } //! returns the spatial size constexpr static Dim_t get_sdim() { return DimS; } //! return a sparse matrix adaptor to the cell Adaptor get_adaptor() override; //! returns the number of degrees of freedom in the cell Dim_t get_nb_dof() const override { return this->size() * ipow(DimS, 2); }; //! return the communicator object const Communicator & get_communicator() const override { return this->projection->get_communicator(); } protected: template Field_t & globalised_field_helper(const std::string & unique_name, int nb_steps_ago); //! make sure that every pixel is assigned to one and only one material void check_material_coverage(); const Ccoord & subdomain_resolutions; //!< the cell's subdomain resolutions const Ccoord & subdomain_locations; //!< the cell's subdomain resolutions const Ccoord & domain_resolutions; //!< the cell's domain resolutions CcoordOps::Pixels pixels; //!< helper to iterate over the pixels const Rcoord & domain_lengths; //!< the cell's lengths Collection_ptr fields; //!< handle for the global fields of the cell StrainField_t & F; //!< ref to strain field StressField_t & P; //!< ref to stress field //! Tangent field might not even be required; so this is an //! optional ref_wrapper instead of a ref optional> K{}; //! container of the materials present in the cell std::vector materials{}; Projection_ptr projection; //!< handle for the projection operator private: }; /** * lightweight resource handle wrapping a `muSpectre::Cell` or * a subclass thereof into `Eigen::EigenBase`, so it can be * interpreted as a sparse matrix by Eigen solvers */ template class CellAdaptor : public Eigen::EigenBase> { public: using Scalar = double; //!< sparse matrix traits using RealScalar = double; //!< sparse matrix traits using StorageIndex = int; //!< sparse matrix traits enum { ColsAtCompileTime = Eigen::Dynamic, MaxColsAtCompileTime = Eigen::Dynamic, RowsAtCompileTime = Eigen::Dynamic, MaxRowsAtCompileTime = Eigen::Dynamic, IsRowMajor = false }; //! constructor explicit CellAdaptor(Cell & cell) : cell{cell} {} //! returns the number of logical rows Eigen::Index rows() const { return this->cell.get_nb_dof(); } //! returns the number of logical columns Eigen::Index cols() const { return this->rows(); } //! implementation of the evaluation template Eigen::Product operator*(const Eigen::MatrixBase & x) const { return Eigen::Product( *this, x.derived()); } Cell & cell; //!< ref to the cell }; } // namespace muSpectre namespace Eigen { namespace internal { //! Implementation of `muSpectre::CellAdaptor` * `Eigen::DenseVector` //! through a specialization of `Eigen::internal::generic_product_impl`: template // GEMV stands for matrix-vector struct generic_product_impl : generic_product_impl_base> { //! undocumented typedef typename Product::Scalar Scalar; //! undocumented template static void scaleAndAddTo(Dest & dst, const CellAdaptor & lhs, const Rhs & rhs, const Scalar & /*alpha*/) { // This method should implement "dst += alpha * lhs * rhs" inplace, // however, for iterative solvers, alpha is always equal to 1, so // let's not bother about it. // Here we could simply call dst.noalias() += lhs.my_matrix() * rhs, dst.noalias() += const_cast(lhs) .cell.evaluate_projected_directional_stiffness(rhs); } }; } // namespace internal } // namespace Eigen #endif // SRC_CELL_CELL_BASE_HH_ diff --git a/src/cell/cell_traits.hh b/src/cell/cell_traits.hh index 7437fae..442e9b0 100644 --- a/src/cell/cell_traits.hh +++ b/src/cell/cell_traits.hh @@ -1,59 +1,58 @@ /** * @file cell_traits.hh * * @author Till Junge * * @date 19 Jan 2018 * * @brief Provides traits for Eigen solvers to be able to use Cells * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #include "common/common.hh" #include #ifndef SRC_CELL_CELL_TRAITS_HH_ #define SRC_CELL_CELL_TRAITS_HH_ namespace muSpectre { - template - class CellAdaptor; + template class CellAdaptor; } // namespace muSpectre namespace Eigen { namespace internal { using Dim_t = muSpectre::Dim_t; //!< universal index type using Real = muSpectre::Real; //!< universal real value type template struct traits> : public Eigen::internal::traits> {}; } // namespace internal } // namespace Eigen #endif // SRC_CELL_CELL_TRAITS_HH_ diff --git a/src/common/field_collection_global.hh b/src/common/field_collection_global.hh index a4290e6..f209413 100644 --- a/src/common/field_collection_global.hh +++ b/src/common/field_collection_global.hh @@ -1,214 +1,216 @@ /** * @file field_collection_global.hh * * @author Till Junge * * @date 05 Nov 2017 * * @brief FieldCollection base-class for global fields * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #ifndef SRC_COMMON_FIELD_COLLECTION_GLOBAL_HH_ #define SRC_COMMON_FIELD_COLLECTION_GLOBAL_HH_ #include "common/ccoord_operations.hh" #include "common/field_collection_base.hh" namespace muSpectre { /** * forward declaration */ template class LocalFieldCollection; /** `GlobalFieldCollection` derives from `FieldCollectionBase` and stores * global fields that live throughout the whole computational domain, i.e. * are defined for each pixel. */ template class GlobalFieldCollection : public FieldCollectionBase> { public: //! for compile time check constexpr static bool Global{true}; using Parent = FieldCollectionBase>; //!< base class //! helpful for functions that fill global fields from local fields using LocalFieldCollection_t = LocalFieldCollection; //! helpful for functions that fill global fields from local fields using GlobalFieldCollection_t = GlobalFieldCollection; using Ccoord = typename Parent::Ccoord; //!< cell coordinates type using Field_p = typename Parent::Field_p; //!< spatial coordinates type //! iterator over all pixels contained it the collection using iterator = typename CcoordOps::Pixels::iterator; //! Default constructor GlobalFieldCollection(); //! Copy constructor GlobalFieldCollection(const GlobalFieldCollection & other) = delete; //! Move constructor GlobalFieldCollection(GlobalFieldCollection && other) = default; //! Destructor virtual ~GlobalFieldCollection() = default; //! Copy assignment operator GlobalFieldCollection & operator=(const GlobalFieldCollection & other) = delete; //! Move assignment operator GlobalFieldCollection & operator=(GlobalFieldCollection && other) = default; /** allocate memory, etc. At this point, the collection is informed aboud the size and shape of the domain (through the sizes parameter). The job of initialise is to make sure that all fields are either of size 0, in which case they need to be allocated, or are of the same size as the product of 'sizes' (if standard strides apply) any field of a different size is wrong. TODO: check whether it makes sense to put a runtime check here **/ inline void initialise(Ccoord sizes, Ccoord locations); //! return subdomain resolutions inline const Ccoord & get_sizes() const; //! return subdomain locations inline const Ccoord & get_locations() const; //! returns the linear index corresponding to cell coordinates template inline size_t get_index(CcoordRef && ccoord) const; //! returns the cell coordinates corresponding to a linear index inline Ccoord get_ccoord(size_t index) const; inline iterator begin() const; //!< returns iterator to first pixel inline iterator end() const; //!< returns iterator past the last pixel //! return spatial dimension (template parameter) static constexpr inline Dim_t spatial_dim() { return DimS; } //! return globalness at compile time static constexpr inline bool is_global() { return Global; } protected: //! number of discretisation cells in each of the DimS spatial directions Ccoord sizes{}; + //! subdomain locations (i.e. coordinates of hind bottom left corner of this + //! subdomain) Ccoord locations{}; CcoordOps::Pixels pixels{}; //!< helper to iterate over the grid private: }; /* ---------------------------------------------------------------------- */ template GlobalFieldCollection::GlobalFieldCollection() : Parent() {} /* ---------------------------------------------------------------------- */ template void GlobalFieldCollection::initialise(Ccoord sizes, Ccoord locations) { if (this->is_initialised) { throw std::runtime_error("double initialisation"); } this->pixels = CcoordOps::Pixels(sizes, locations); this->size_ = CcoordOps::get_size(sizes); this->sizes = sizes; this->locations = locations; std::for_each( std::begin(this->fields), std::end(this->fields), [this](auto && item) { auto && field = *item.second; const auto field_size = field.size(); if (field_size == 0) { field.resize(this->size()); } else if (field_size != this->size()) { std::stringstream err_stream; err_stream << "Field '" << field.get_name() << "' contains " << field_size << " entries, but the field collection " << "has " << this->size() << " pixels"; throw FieldCollectionError(err_stream.str()); } }); this->is_initialised = true; } //----------------------------------------------------------------------------// //! return subdomain resolutions template const typename GlobalFieldCollection::Ccoord & GlobalFieldCollection::get_sizes() const { return this->sizes; } //----------------------------------------------------------------------------// //! return subdomain locations template const typename GlobalFieldCollection::Ccoord & GlobalFieldCollection::get_locations() const { return this->locations; } //----------------------------------------------------------------------------// //! returns the cell coordinates corresponding to a linear index template typename GlobalFieldCollection::Ccoord GlobalFieldCollection::get_ccoord(size_t index) const { return CcoordOps::get_ccoord(this->get_sizes(), this->get_locations(), std::move(index)); } /* ---------------------------------------------------------------------- */ template typename GlobalFieldCollection::iterator GlobalFieldCollection::begin() const { return this->pixels.begin(); } /* ---------------------------------------------------------------------- */ template typename GlobalFieldCollection::iterator GlobalFieldCollection::end() const { return this->pixels.end(); } //-------------------------------------------------------------------------// //! returns the linear index corresponding to cell coordinates template template size_t GlobalFieldCollection::get_index(CcoordRef && ccoord) const { static_assert( std::is_same>>::value, "can only be called with values or references of Ccoord"); return CcoordOps::get_index(this->get_sizes(), this->get_locations(), std::forward(ccoord)); } } // namespace muSpectre #endif // SRC_COMMON_FIELD_COLLECTION_GLOBAL_HH_ diff --git a/src/materials/material_evaluator.hh b/src/materials/material_evaluator.hh new file mode 100644 index 0000000..32866ce --- /dev/null +++ b/src/materials/material_evaluator.hh @@ -0,0 +1,267 @@ +/** + * @file material_evaluator.hh + * + * @author Till Junge + * + * @date 12 Dec 2018 + * + * @brief Helper to evaluate material laws + * + * Copyright © 2018 Till Junge + * + * µSpectre is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License as + * published by the Free Software Foundation, either version 3, or (at + * your option) any later version. + * + * µSpectre 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 Lesser General Public License + * along with µSpectre; see the file COPYING. If not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + * + * Additional permission under GNU GPL version 3 section 7 + * + * If you modify this Program, or any covered work, by linking or combining it + * with proprietary FFT implementations or numerical libraries, containing parts + * covered by the terms of those libraries' licenses, the licensors of this + * Program grant you additional permission to convey the resulting work. + */ + +#ifndef SRC_MATERIALS_MATERIAL_EVALUATOR_HH_ +#define SRC_MATERIALS_MATERIAL_EVALUATOR_HH_ + +#include "common/common.hh" +#include "common/T4_map_proxy.hh" +#include "common/ccoord_operations.hh" +#include "materials/materials_toolbox.hh" +#include "common/field.hh" + +#include +#include +#include + +namespace muSpectre { + + /** + * forward declaration to avoid incluning material_base.hh + */ + template + class MaterialBase; + + template + class MaterialEvaluator { + public: + using T2_t = Eigen::Matrix; + using T4_t = T4Mat; + + using T2_map = Eigen::Map; + using T4_map = T4MatMap; + + using T2_const_map = Eigen::Map; + using T4_const_map = T4MatMap; + + using FieldColl_t = GlobalFieldCollection; + using T2Field_t = TensorField; + using T4Field_t = TensorField; + + //! Default constructor + MaterialEvaluator() = delete; + + /** + * constructor with a shared pointer to a Material + */ + explicit MaterialEvaluator( + std::shared_ptr> material) + : material{material}, + collection{std::make_unique()}, + strain{make_field("gradient", *this->collection)}, + stress{make_field("stress", *this->collection)}, + tangent{make_field("tangent", *this->collection)} { + this->collection->initialise(CcoordOps::get_cube(1), {0}); + } + + //! Copy constructor + MaterialEvaluator(const MaterialEvaluator & other) = delete; + + //! Move constructor + MaterialEvaluator(MaterialEvaluator && other) = default; + + //! Destructor + virtual ~MaterialEvaluator() = default; + + //! Copy assignment operator + MaterialEvaluator & operator=(const MaterialEvaluator & other) = delete; + + //! Move assignment operator + MaterialEvaluator & operator=(MaterialEvaluator && other) = default; + + /** + * for materials with state variables. See `muSpectre::MaterialBase` for + * details + */ + void save_history_variables() { this->material->save_history_variables(); } + + /** + * Evaluates the underlying materials constitutive law and returns the + * stress P or σ as a function of the placement gradient F or small strain + * tensor ε depending on the formulation + * (`muSpectre::Formulation::small_strain` for σ(ε), + * `muSpectre::Formulation::finite_strain` for P(F)) + */ + inline T2_const_map evaluate_stress(const Eigen::Ref & grad, + const Formulation & form); + + /** + * Evaluates the underlying materials constitutive law and returns the the + * stress P or σ and the tangent moduli K as a function of the placement + * gradient F or small strain tensor ε depending on the formulation + * (`muSpectre::Formulation::small_strain` for σ(ε), + * `muSpectre::Formulation::finite_strain` for P(F)) + */ + inline std::tuple + evaluate_stress_tangent(const Eigen::Ref & grad, + const Formulation & form); + + /** + * estimate the tangent using finite difference + */ + inline T4_t + estimate_tangent(const Eigen::Ref & grad, + const Formulation & form, const Real step, + const FiniteDiff diff_type = FiniteDiff::centred); + + protected: + /** + * throws a runtime error if the material's per-pixel data has not been set. + */ + void check_init() const; + + std::shared_ptr> material; + std::unique_ptr collection; + T2Field_t & strain; + T2Field_t & stress; + T4Field_t & tangent; + }; + + /* ---------------------------------------------------------------------- */ + template + auto MaterialEvaluator::evaluate_stress( + const Eigen::Ref & grad, const Formulation & form) + -> T2_const_map { + this->check_init(); + this->strain.get_map()[0] = grad; + this->material->compute_stresses(this->strain, this->stress, form); + return T2_const_map(this->stress.get_map()[0].data()); + } + + /* ---------------------------------------------------------------------- */ + template + auto MaterialEvaluator::evaluate_stress_tangent( + const Eigen::Ref & grad, const Formulation & form) + -> std::tuple { + this->check_init(); + this->strain.get_map()[0] = grad; + this->material->compute_stresses_tangent(this->strain, this->stress, + this->tangent, form); + return std::make_tuple(T2_const_map(this->stress.get_map()[0].data()), + T4_const_map(this->tangent.get_map()[0].data())); + } + + template + void MaterialEvaluator::check_init() const { + const auto & size{this->material->size()}; + if (size < 1) { + throw std::runtime_error( + "Not initialised! You have to call `add_pixel(...)` on your material " + "exactly one time before you can evaluate it."); + } else if (size > 1) { + std::stringstream error{}; + error << "The material to be evaluated should have exactly one pixel " + "added. You've added " + << size << " pixels."; + throw std::runtime_error(error.str()); + } + } + + /* ---------------------------------------------------------------------- */ + template + auto MaterialEvaluator::estimate_tangent( + const Eigen::Ref & grad, const Formulation & form, + const Real delta, const FiniteDiff diff_type) -> T4_t { + using T2_vec = Eigen::Map>; + + T4_t tangent{T4_t::Zero()}; + + const T2_t stress{this->evaluate_stress(grad, form)}; + + auto fun{[&form, this](const auto & grad_) { + return this->evaluate_stress(grad_, form); + }}; + + auto symmetrise{ + [](auto & eps) { eps = .5 * (eps + eps.transpose().eval()); }}; + + static_assert(Int(decltype(tangent.col(0))::SizeAtCompileTime) == + Int(T2_t::SizeAtCompileTime), + "wrong column size"); + + for (Dim_t i{}; i < DimM * DimM; ++i) { + T2_t strain2{grad}; + T2_vec strain2_vec{strain2.data()}; + switch (diff_type) { + case FiniteDiff::forward: { + strain2_vec(i) += delta; + if (form == Formulation::small_strain) { + symmetrise(strain2); + } + + T2_t del_f_del{(fun(strain2) - stress) / delta}; + + tangent.col(i) = T2_vec(del_f_del.data()); + break; + } + case FiniteDiff::backward: { + strain2_vec(i) -= delta; + if (form == Formulation::small_strain) { + symmetrise(strain2); + } + T2_t del_f_del{(stress - fun(strain2)) / delta}; + + tangent.col(i) = T2_vec(del_f_del.data()); + break; + } + case FiniteDiff::centred: { + T2_t strain1{grad}; + T2_vec strain1_vec{strain1.data()}; + strain1_vec(i) += delta; + strain2_vec(i) -= delta; + if (form == Formulation::small_strain) { + symmetrise(strain1); + symmetrise(strain2); + } + T2_t del_f_del{(fun(strain1).eval() - fun(strain2).eval()) / + (2 * delta)}; + + tangent.col(i) = T2_vec(del_f_del.data()); + break; + } + default: { + throw std::runtime_error("Unknown FiniteDiff value"); + break; + } + } + static_assert(Int(decltype(tangent.col(i))::SizeAtCompileTime) == + Int(T2_t::SizeAtCompileTime), + "wrong column size"); + } + return tangent; + } + +} // namespace muSpectre + +#endif // SRC_MATERIALS_MATERIAL_EVALUATOR_HH_ diff --git a/src/materials/material_linear_elastic1.hh b/src/materials/material_linear_elastic1.hh index 50642e7..4a06a20 100644 --- a/src/materials/material_linear_elastic1.hh +++ b/src/materials/material_linear_elastic1.hh @@ -1,195 +1,188 @@ /** * @file material_linear_elastic1.hh * * @author Till Junge * * @date 13 Nov 2017 * * @brief Implementation for linear elastic reference material like in de Geus * 2017. This follows the simplest and likely not most efficient * implementation (with exception of the Python law) * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #ifndef SRC_MATERIALS_MATERIAL_LINEAR_ELASTIC1_HH_ #define SRC_MATERIALS_MATERIAL_LINEAR_ELASTIC1_HH_ #include "common/common.hh" #include "materials/stress_transformations_PK2.hh" #include "materials/material_muSpectre_base.hh" #include "materials/materials_toolbox.hh" namespace muSpectre { template class MaterialLinearElastic1; /** * traits for objective linear elasticity */ template struct MaterialMuSpectre_traits> { using Parent = MaterialMuSpectre_traits; //!< base for elasticity //! global field collection using GFieldCollection_t = typename MaterialBase::GFieldCollection_t; //! expected map type for strain fields using StrainMap_t = MatrixFieldMap; //! expected map type for stress fields using StressMap_t = MatrixFieldMap; //! expected map type for tangent stiffness fields using TangentMap_t = T4MatrixFieldMap; //! declare what type of strain measure your law takes as input constexpr static auto strain_measure{StrainMeasure::GreenLagrange}; //! declare what type of stress measure your law yields as output constexpr static auto stress_measure{StressMeasure::PK2}; //! elasticity without internal variables using InternalVariables = std::tuple<>; }; //! DimS spatial dimension (dimension of problem //! DimM material_dimension (dimension of constitutive law) /** * implements objective linear elasticity */ template class MaterialLinearElastic1 : public MaterialMuSpectre, DimS, DimM> { public: //! base class using Parent = MaterialMuSpectre; /** * type used to determine whether the * `muSpectre::MaterialMuSpectre::iterable_proxy` evaluate only * stresses or also tangent stiffnesses */ using NeedTangent = typename Parent::NeedTangent; //! global field collection using Stiffness_t = Eigen::TensorFixedSize>; //! traits of this material using traits = MaterialMuSpectre_traits; //! this law does not have any internal variables using InternalVariables = typename traits::InternalVariables; //! Hooke's law implementation using Hooke = typename MatTB::Hooke; //! Default constructor MaterialLinearElastic1() = delete; //! Copy constructor MaterialLinearElastic1(const MaterialLinearElastic1 & other) = delete; //! Construct by name, Young's modulus and Poisson's ratio MaterialLinearElastic1(std::string name, Real young, Real poisson); //! Move constructor MaterialLinearElastic1(MaterialLinearElastic1 && other) = delete; //! Destructor virtual ~MaterialLinearElastic1() = default; //! Copy assignment operator MaterialLinearElastic1 & operator=(const MaterialLinearElastic1 & other) = delete; //! Move assignment operator MaterialLinearElastic1 & operator=(MaterialLinearElastic1 && other) = delete; /** * evaluates second Piola-Kirchhoff stress given the Green-Lagrange * strain (or Cauchy stress if called with a small strain tensor) */ template inline decltype(auto) evaluate_stress(s_t && E); /** * evaluates both second Piola-Kirchhoff stress and stiffness given * the Green-Lagrange strain (or Cauchy stress and stiffness if * called with a small strain tensor) */ template inline decltype(auto) evaluate_stress_tangent(s_t && E); /** * return the empty internals tuple */ InternalVariables & get_internals() { return this->internal_variables; } protected: const Real young; //!< Young's modulus const Real poisson; //!< Poisson's ratio const Real lambda; //!< first Lamé constant const Real mu; //!< second Lamé constant (shear modulus) const Stiffness_t C; //!< stiffness tensor //! empty tuple InternalVariables internal_variables{}; private: }; /* ---------------------------------------------------------------------- */ template template auto MaterialLinearElastic1::evaluate_stress(s_t && E) -> decltype(auto) { return Hooke::evaluate_stress(this->lambda, this->mu, std::move(E)); } /* ---------------------------------------------------------------------- */ template template auto MaterialLinearElastic1::evaluate_stress_tangent(s_t && E) -> decltype(auto) { using Tangent_t = typename traits::TangentMap_t::reference; - // using mat = Eigen::Matrix; - // mat ecopy{E}; - // std::cout << "E" << std::endl << ecopy << std::endl; - // std::cout << "P1" << std::endl << mat{ - // std::get<0>(Hooke::evaluate_stress(this->lambda, this->mu, - // Tangent_t(const_cast(this->C.data())), - // std::move(E)))} << std::endl; return Hooke::evaluate_stress( this->lambda, this->mu, Tangent_t(const_cast(this->C.data())), std::move(E)); } } // namespace muSpectre #endif // SRC_MATERIALS_MATERIAL_LINEAR_ELASTIC1_HH_ diff --git a/src/materials/material_linear_elastic2.hh b/src/materials/material_linear_elastic2.hh index 873c33b..b9a4e6f 100644 --- a/src/materials/material_linear_elastic2.hh +++ b/src/materials/material_linear_elastic2.hh @@ -1,201 +1,201 @@ /** * @file material_linear_elastic2.hh * * @author Till Junge * * @date 03 Feb 2018 * * @brief linear elastic material with imposed eigenstrain and its * type traits. Uses the MaterialMuSpectre facilities to keep it * simple * * Copyright © 2018 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #ifndef SRC_MATERIALS_MATERIAL_LINEAR_ELASTIC2_HH_ #define SRC_MATERIALS_MATERIAL_LINEAR_ELASTIC2_HH_ #include "materials/material_linear_elastic1.hh" #include "common/field.hh" #include namespace muSpectre { template class MaterialLinearElastic2; /** * traits for objective linear elasticity with eigenstrain */ template struct MaterialMuSpectre_traits> { //! global field collection using GFieldCollection_t = typename MaterialBase::GFieldCollection_t; //! expected map type for strain fields using StrainMap_t = MatrixFieldMap; //! expected map type for stress fields using StressMap_t = MatrixFieldMap; //! expected map type for tangent stiffness fields using TangentMap_t = T4MatrixFieldMap; //! declare what type of strain measure your law takes as input constexpr static auto strain_measure{StrainMeasure::GreenLagrange}; //! declare what type of stress measure your law yields as output constexpr static auto stress_measure{StressMeasure::PK2}; //! local field_collections used for internals using LFieldColl_t = LocalFieldCollection; //! local strain type using LStrainMap_t = MatrixFieldMap; //! elasticity with eigenstrain using InternalVariables = std::tuple; }; /** * implements objective linear elasticity with an eigenstrain per pixel */ template class MaterialLinearElastic2 : public MaterialMuSpectre, DimS, DimM> { public: //! base class using Parent = MaterialMuSpectre; //! type for stiffness tensor construction using Stiffness_t = Eigen::TensorFixedSize>; //! traits of this material using traits = MaterialMuSpectre_traits; //! Type of container used for storing eigenstrain using InternalVariables = typename traits::InternalVariables; //! Hooke's law implementation using Hooke = typename MatTB::Hooke; //! reference to any type that casts to a matrix - using StrainTensor = Eigen::Ref>; + using StrainTensor = Eigen::Ref>; //! Default constructor MaterialLinearElastic2() = delete; //! Construct by name, Young's modulus and Poisson's ratio MaterialLinearElastic2(std::string name, Real young, Real poisson); //! Copy constructor MaterialLinearElastic2(const MaterialLinearElastic2 & other) = delete; //! Move constructor MaterialLinearElastic2(MaterialLinearElastic2 && other) = delete; //! Destructor virtual ~MaterialLinearElastic2() = default; //! Copy assignment operator MaterialLinearElastic2 & operator=(const MaterialLinearElastic2 & other) = delete; //! Move assignment operator MaterialLinearElastic2 & operator=(MaterialLinearElastic2 && other) = delete; /** * evaluates second Piola-Kirchhoff stress given the Green-Lagrange * strain (or Cauchy stress if called with a small strain tensor) */ template inline decltype(auto) evaluate_stress(s_t && E, eigen_s_t && E_eig); /** * evaluates both second Piola-Kirchhoff stress and stiffness given * the Green-Lagrange strain (or Cauchy stress and stiffness if * called with a small strain tensor) */ template inline decltype(auto) evaluate_stress_tangent(s_t && E, eigen_s_t && E_eig); /** * return the internals tuple */ InternalVariables & get_internals() { return this->internal_variables; } /** * overload add_pixel to write into eigenstrain */ void add_pixel(const Ccoord_t & pixel) final; /** * overload add_pixel to write into eigenstrain */ void add_pixel(const Ccoord_t & pixel, const StrainTensor & E_eig); protected: //! linear material without eigenstrain used to compute response MaterialLinearElastic1 material; //! storage for eigenstrain using Field_t = TensorField, Real, secondOrder, DimM>; Field_t & eigen_field; //!< field holding the eigen strain per pixel //! tuple for iterable eigen_field InternalVariables internal_variables; }; /* ---------------------------------------------------------------------- */ template template auto MaterialLinearElastic2::evaluate_stress(s_t && E, eigen_s_t && E_eig) -> decltype(auto) { return this->material.evaluate_stress(E - E_eig); } /* ---------------------------------------------------------------------- */ template template auto MaterialLinearElastic2::evaluate_stress_tangent( s_t && E, eigen_s_t && E_eig) -> decltype(auto) { // using mat = Eigen::Matrix; // mat ecopy{E}; // mat eig_copy{E_eig}; // mat ediff{ecopy-eig_copy}; // std::cout << "eidff - (E-E_eig)" << std::endl << ediff-(E-E_eig) << // std::endl; std::cout << "P1 " << std::endl << // mat{std::get<0>(this->material.evaluate_stress_tangent(E-E_eig))} << // "" << std::endl; std::cout << "P2" << std::endl << // mat{std::get<0>(this->material.evaluate_stress_tangent(std::move(ediff)))} // << std::endl; return this->material.evaluate_stress_tangent(E - E_eig); } } // namespace muSpectre #endif // SRC_MATERIALS_MATERIAL_LINEAR_ELASTIC2_HH_ diff --git a/src/materials/material_muSpectre_base.hh b/src/materials/material_muSpectre_base.hh index 14620a5..491c3bb 100644 --- a/src/materials/material_muSpectre_base.hh +++ b/src/materials/material_muSpectre_base.hh @@ -1,699 +1,721 @@ /** * @file material_muSpectre_base.hh * * @author Till Junge * * @date 25 Oct 2017 * * @brief Base class for materials written for µSpectre specifically. These * can take full advantage of the configuration-change utilities of * µSpectre. The user can inherit from them to define new constitutive * laws and is merely required to provide the methods for computing the * second Piola-Kirchhoff stress and Tangent. This class uses the * "curiously recurring template parameter" to avoid virtual calls. * * Copyright © 2017 Till Junge * * µSpectre is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public License as * published by the Free Software Foundation, either version 3, or (at * your option) any later version. * * µSpectre 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 Lesser General Public License * along with µSpectre; see the file COPYING. If not, write to the * Free Software Foundation, Inc., 59 Temple Place - Suite 330, * * Boston, MA 02111-1307, USA. * * Additional permission under GNU GPL version 3 section 7 * * If you modify this Program, or any covered work, by linking or combining it * with proprietary FFT implementations or numerical libraries, containing parts * covered by the terms of those libraries' licenses, the licensors of this * Program grant you additional permission to convey the resulting work. */ #ifndef SRC_MATERIALS_MATERIAL_MUSPECTRE_BASE_HH_ #define SRC_MATERIALS_MATERIAL_MUSPECTRE_BASE_HH_ #include "common/common.hh" #include "materials/material_base.hh" #include "materials/materials_toolbox.hh" +#include "materials/material_evaluator.hh" #include "common/field_collection.hh" #include "common/field.hh" #include "common//utilities.hh" #include #include #include #include namespace muSpectre { // Forward declaration for factory function template class CellBase; /** * material traits are used by `muSpectre::MaterialMuSpectre` to * break the circular dependence created by the curiously recurring * template parameter. These traits must define * - these `muSpectre::FieldMap`s: * - `StrainMap_t`: typically a `muSpectre::MatrixFieldMap` for a * constant second-order `muSpectre::TensorField` * - `StressMap_t`: typically a `muSpectre::MatrixFieldMap` for a * writable secord-order `muSpectre::TensorField` * - `TangentMap_t`: typically a `muSpectre::T4MatrixFieldMap` for a * writable fourth-order `muSpectre::TensorField` * - `strain_measure`: the expected strain type (will be replaced by the * small-strain tensor ε * `muspectre::StrainMeasure::Infinitesimal` in small * strain computations) * - `stress_measure`: the measure of the returned stress. Is used by * `muspectre::MaterialMuSpectre` to transform it into * Cauchy stress (`muspectre::StressMeasure::Cauchy`) in * small-strain computations and into first * Piola-Kirchhoff stress `muspectre::StressMeasure::PK1` * in finite-strain computations * - `InternalVariables`: a tuple of `muSpectre::FieldMap`s containing * internal variables */ template struct MaterialMuSpectre_traits {}; template class MaterialMuSpectre; /** * Base class for most convenient implementation of materials */ template class MaterialMuSpectre : public MaterialBase { public: /** * type used to determine whether the * `muSpectre::MaterialMuSpectre::iterable_proxy` evaluate only * stresses or also tangent stiffnesses */ using NeedTangent = MatTB::NeedTangent; using Parent = MaterialBase; //!< base class //! global field collection using GFieldCollection_t = typename Parent::GFieldCollection_t; //! expected type for stress fields using StressField_t = typename Parent::StressField_t; //! expected type for strain fields using StrainField_t = typename Parent::StrainField_t; //! expected type for tangent stiffness fields using TangentField_t = typename Parent::TangentField_t; //! traits for the CRTP subclass using traits = MaterialMuSpectre_traits; //! Default constructor MaterialMuSpectre() = delete; //! Construct by name explicit MaterialMuSpectre(std::string name); //! Copy constructor MaterialMuSpectre(const MaterialMuSpectre & other) = delete; //! Move constructor MaterialMuSpectre(MaterialMuSpectre && other) = delete; //! Destructor virtual ~MaterialMuSpectre() = default; //! Factory template static Material & make(CellBase & cell, ConstructorArgs &&... args); + /** Factory + * takes all arguments after the name of the underlying + * Material's constructor. E.g., if the underlying material is a + * `muSpectre::MaterialLinearElastic1`, these would be Young's + * modulus and Poisson's ratio. + */ + template + static std::tuple, MaterialEvaluator> + make_evaluator(ConstructorArgs &&... args); //! Copy assignment operator MaterialMuSpectre & operator=(const MaterialMuSpectre & other) = delete; //! Move assignment operator MaterialMuSpectre & operator=(MaterialMuSpectre && other) = delete; //! allocate memory, etc void initialise() override; using Parent::compute_stresses; using Parent::compute_stresses_tangent; //! computes stress void compute_stresses(const StrainField_t & F, StressField_t & P, Formulation form) final; //! computes stress and tangent modulus void compute_stresses_tangent(const StrainField_t & F, StressField_t & P, TangentField_t & K, Formulation form) final; protected: //! computes stress with the formulation available at compile time //! __attribute__ required by g++-6 and g++-7 because of this bug: //! https://gcc.gnu.org/bugzilla/show_bug.cgi?id=80947 template inline void compute_stresses_worker(const StrainField_t & F, StressField_t & P) __attribute__((visibility("default"))); //! computes stress with the formulation available at compile time //! __attribute__ required by g++-6 and g++-7 because of this bug: //! https://gcc.gnu.org/bugzilla/show_bug.cgi?id=80947 template inline void compute_stresses_worker(const StrainField_t & F, StressField_t & P, TangentField_t & K) __attribute__((visibility("default"))); //! this iterable class is a default for simple laws that just take a strain //! the iterable is just a templated wrapper to provide a range to iterate //! over that does or does not include tangent moduli template class iterable_proxy; /** * inheriting classes with internal variables need to overload this function */ typename traits::InternalVariables & get_internals() { return static_cast(*this).get_internals(); } bool is_initialised{false}; //!< to handle double initialisation right private: }; /* ---------------------------------------------------------------------- */ template MaterialMuSpectre::MaterialMuSpectre(std::string name) : Parent(name) { using stress_compatible = typename traits::StressMap_t::template is_compatible; using strain_compatible = typename traits::StrainMap_t::template is_compatible; using tangent_compatible = typename traits::TangentMap_t::template is_compatible; static_assert((stress_compatible::value && stress_compatible::explain()), "The material's declared stress map is not compatible " "with the stress field. More info in previously shown " "assert."); static_assert((strain_compatible::value && strain_compatible::explain()), "The material's declared strain map is not compatible " "with the strain field. More info in previously shown " "assert."); static_assert((tangent_compatible::value && tangent_compatible::explain()), "The material's declared tangent map is not compatible " "with the tangent field. More info in previously shown " "assert."); } /* ---------------------------------------------------------------------- */ template template Material & MaterialMuSpectre::make(CellBase & cell, ConstructorArgs &&... args) { auto mat = std::make_unique(args...); auto & mat_ref = *mat; cell.add_material(std::move(mat)); return mat_ref; } + /* ---------------------------------------------------------------------- */ + template + template + std::tuple, MaterialEvaluator> + MaterialMuSpectre::make_evaluator( + ConstructorArgs &&... args) { + auto mat = std::make_shared("name", args...); + using Ret_t = + std::tuple, MaterialEvaluator>; + return Ret_t(mat, MaterialEvaluator{mat}); + } + /* ---------------------------------------------------------------------- */ template void MaterialMuSpectre::initialise() { if (!this->is_initialised) { this->internal_fields.initialise(); this->is_initialised = true; } } /* ---------------------------------------------------------------------- */ template void MaterialMuSpectre::compute_stresses( const StrainField_t & F, StressField_t & P, Formulation form) { switch (form) { case Formulation::finite_strain: { this->template compute_stresses_worker(F, P); break; } case Formulation::small_strain: { this->template compute_stresses_worker(F, P); break; } default: throw std::runtime_error("Unknown formulation"); break; } } /* ---------------------------------------------------------------------- */ template void MaterialMuSpectre::compute_stresses_tangent( const StrainField_t & F, StressField_t & P, TangentField_t & K, Formulation form) { switch (form) { case Formulation::finite_strain: { this->template compute_stresses_worker(F, P, K); break; } case Formulation::small_strain: { this->template compute_stresses_worker(F, P, K); break; } default: throw std::runtime_error("Unknown formulation"); break; } } /* ---------------------------------------------------------------------- */ template template void MaterialMuSpectre::compute_stresses_worker( const StrainField_t & F, StressField_t & P, TangentField_t & K) { /* These lambdas are executed for every integration point. F contains the transformation gradient for finite strain calculations and the infinitesimal strain tensor in small strain problems The internal_variables tuple contains whatever internal variables Material declared (e.g., eigenstrain, strain rate, etc.) */ using Strains_t = std::tuple; using Stresses_t = std::tuple; auto constitutive_law_small_strain = [this](Strains_t Strains, Stresses_t Stresses, auto && internal_variables) { constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; constexpr StrainMeasure expected_strain_m{ get_formulation_strain_type(Form, traits::strain_measure)}; auto & this_mat = static_cast(*this); // Transformation gradient is first in the strains tuple auto & F = std::get<0>(Strains); auto && strain = MatTB::convert_strain(F); // return value contains a tuple of rvalue_refs to both stress and tangent // moduli Stresses = apply( [&strain, &this_mat](auto &&... internals) { return this_mat.evaluate_stress_tangent(std::move(strain), internals...); }, internal_variables); }; auto constitutive_law_finite_strain = [this](Strains_t Strains, Stresses_t Stresses, auto && internal_variables) { constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; constexpr StrainMeasure expected_strain_m{ get_formulation_strain_type(Form, traits::strain_measure)}; auto & this_mat = static_cast(*this); // Transformation gradient is first in the strains tuple auto & grad = std::get<0>(Strains); auto && strain = MatTB::convert_strain(grad); // TODO(junge): Figure this out: I can't std::move(internals...), // because if there are no internals, compilation fails with "no // matching function for call to ‘move()’'. These are tuples of // lvalue references, so it shouldn't be too bad, but still // irksome. // return value contains a tuple of rvalue_refs to both stress // and tangent moduli auto stress_tgt = apply( [&strain, &this_mat](auto &&... internals) { return this_mat.evaluate_stress_tangent(std::move(strain), internals...); }, internal_variables); auto & stress = std::get<0>(stress_tgt); auto & tangent = std::get<1>(stress_tgt); Stresses = MatTB::PK1_stress( std::move(grad), std::move(stress), std::move(tangent)); }; iterable_proxy fields{*this, F, P, K}; for (auto && arglist : fields) { /** * arglist is a tuple of three tuples containing only Lvalue * references (see value_tye in the class definition of * iterable_proxy::iterator). Tuples contain strains, stresses * and internal variables, respectively, */ // auto && stress_tgt = std::get<0>(tuples); // auto && inputs = std::get<1>(tuples);TODO:clean this static_assert(std::is_same( std::get<0>(arglist)))>>::value, "Type mismatch for strain reference, check iterator " "value_type"); static_assert(std::is_same( std::get<1>(arglist)))>>::value, "Type mismatch for stress reference, check iterator" "value_type"); static_assert(std::is_same( std::get<1>(arglist)))>>::value, "Type mismatch for tangent reference, check iterator" "value_type"); switch (Form) { case Formulation::small_strain: { apply(constitutive_law_small_strain, std::move(arglist)); break; } case Formulation::finite_strain: { apply(constitutive_law_finite_strain, std::move(arglist)); break; } } } } /* ---------------------------------------------------------------------- */ template template void MaterialMuSpectre::compute_stresses_worker( const StrainField_t & F, StressField_t & P) { /* These lambdas are executed for every integration point. F contains the transformation gradient for finite strain calculations and the infinitesimal strain tensor in small strain problems The internal_variables tuple contains whatever internal variables Material declared (e.g., eigenstrain, strain rate, etc.) */ using Strains_t = std::tuple; using Stresses_t = std::tuple; auto constitutive_law_small_strain = [this](Strains_t Strains, Stresses_t Stresses, auto && internal_variables) { constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; constexpr StrainMeasure expected_strain_m{ get_formulation_strain_type(Form, traits::strain_measure)}; auto & this_mat = static_cast(*this); // Transformation gradient is first in the strains tuple auto & F = std::get<0>(Strains); auto && strain = MatTB::convert_strain(F); // return value contains a tuple of rvalue_refs to both stress and tangent // moduli auto & sigma = std::get<0>(Stresses); sigma = apply( [&strain, &this_mat](auto &&... internals) { return this_mat.evaluate_stress(std::move(strain), internals...); }, internal_variables); }; auto constitutive_law_finite_strain = [this](Strains_t Strains, Stresses_t && Stresses, auto && internal_variables) { constexpr StrainMeasure stored_strain_m{get_stored_strain_type(Form)}; constexpr StrainMeasure expected_strain_m{ get_formulation_strain_type(Form, traits::strain_measure)}; auto & this_mat = static_cast(*this); // Transformation gradient is first in the strains tuple auto & F = std::get<0>(Strains); auto && strain = MatTB::convert_strain(F); // TODO(junge): Figure this out: I can't std::move(internals...), // because if there are no internals, compilation fails with "no // matching function for call to ‘move()’'. These are tuples of // lvalue references, so it shouldn't be too bad, but still // irksome. // return value contains a tuple of rvalue_refs to both stress // and tangent moduli auto && stress = apply( [&strain, &this_mat](auto &&... internals) { return this_mat.evaluate_stress(std::move(strain), internals...); }, internal_variables); auto & P = get<0>(Stresses); P = MatTB::PK1_stress( F, stress); }; iterable_proxy fields{*this, F, P}; for (auto && arglist : fields) { /** * arglist is a tuple of three tuples containing only Lvalue * references (see value_tye in the class definition of * iterable_proxy::iterator). Tuples contain strains, stresses * and internal variables, respectively, */ // auto && stress_tgt = std::get<0>(tuples); // auto && inputs = std::get<1>(tuples);TODO:clean this static_assert(std::is_same( std::get<0>(arglist)))>>::value, "Type mismatch for strain reference, check iterator " "value_type"); static_assert(std::is_same( std::get<1>(arglist)))>>::value, "Type mismatch for stress reference, check iterator" "value_type"); switch (Form) { case Formulation::small_strain: { apply(constitutive_law_small_strain, std::move(arglist)); break; } case Formulation::finite_strain: { apply(constitutive_law_finite_strain, std::move(arglist)); break; } } } } /* ---------------------------------------------------------------------- */ //! this iterator class is a default for simple laws that just take a strain template template class MaterialMuSpectre::iterable_proxy { public: //! Default constructor iterable_proxy() = delete; /** * type used to determine whether the * `muSpectre::MaterialMuSpectre::iterable_proxy` evaluate only * stresses or also tangent stiffnesses */ using NeedTangent = typename MaterialMuSpectre::NeedTangent; /** Iterator uses the material's internal variables field collection to iterate selectively over the global fields (such as the transformation gradient F and first Piola-Kirchhoff stress P. **/ template iterable_proxy(MaterialMuSpectre & mat, const StrainField_t & F, StressField_t & P, std::enable_if_t & K) : material{mat}, strain_field{F}, stress_tup{P, K}, internals{material.get_internals()} {}; /** Iterator uses the material's internal variables field collection to iterate selectively over the global fields (such as the transformation gradient F and first Piola-Kirchhoff stress P. **/ template iterable_proxy(MaterialMuSpectre & mat, const StrainField_t & F, std::enable_if_t & P) : material{mat}, strain_field{F}, stress_tup{P}, internals{material.get_internals()} {}; //! Expected type for strain fields using StrainMap_t = typename traits::StrainMap_t; //! Expected type for stress fields using StressMap_t = typename traits::StressMap_t; //! Expected type for tangent stiffness fields using TangentMap_t = typename traits::TangentMap_t; //! expected type for strain values using Strain_t = typename traits::StrainMap_t::reference; //! expected type for stress values using Stress_t = typename traits::StressMap_t::reference; //! expected type for tangent stiffness values using Tangent_t = typename traits::TangentMap_t::reference; //! tuple of intervnal variables, depends on the material using InternalVariables = typename traits::InternalVariables; //! tuple containing a stress and possibly a tangent stiffness field using StressFieldTup = std::conditional_t<(NeedTgt == NeedTangent::yes), std::tuple, std::tuple>; //! tuple containing a stress and possibly a tangent stiffness field map using StressMapTup = std::conditional_t<(NeedTgt == NeedTangent::yes), std::tuple, std::tuple>; //! tuple containing a stress and possibly a tangent stiffness value ref using Stress_tTup = std::conditional_t<(NeedTgt == NeedTangent::yes), std::tuple, std::tuple>; //! Copy constructor iterable_proxy(const iterable_proxy & other) = default; //! Move constructor iterable_proxy(iterable_proxy && other) = default; //! Destructor virtual ~iterable_proxy() = default; //! Copy assignment operator iterable_proxy & operator=(const iterable_proxy & other) = default; //! Move assignment operator iterable_proxy & operator=(iterable_proxy && other) = default; /** * dereferences into a tuple containing strains, and internal * variables, as well as maps to the stress and potentially * stiffness maps where to write the response of a pixel */ class iterator { public: //! type to refer to internal variables owned by a CRTP material using InternalReferences = MatTB::ReferenceTuple_t; //! return type to be unpacked per pixel my the constitutive law using value_type = std::tuple, Stress_tTup, InternalReferences>; using iterator_category = std::forward_iterator_tag; //!< stl conformance //! Default constructor iterator() = delete; /** Iterator uses the material's internal variables field collection to iterate selectively over the global fields (such as the transformation gradient F and first Piola-Kirchhoff stress P. **/ explicit iterator(const iterable_proxy & it, bool begin = true) : it{it}, strain_map{it.strain_field}, stress_map(it.stress_tup), index{begin ? 0 : it.material.internal_fields.size()} {} //! Copy constructor iterator(const iterator & other) = default; //! Move constructor iterator(iterator && other) = default; //! Destructor virtual ~iterator() = default; //! Copy assignment operator iterator & operator=(const iterator & other) = default; //! Move assignment operator iterator & operator=(iterator && other) = default; //! pre-increment inline iterator & operator++(); //! dereference inline value_type operator*(); //! inequality inline bool operator!=(const iterator & other) const; protected: const iterable_proxy & it; //!< ref to the proxy StrainMap_t strain_map; //!< map onto the global strain field //! map onto the global stress field and possibly tangent stiffness StressMapTup stress_map; size_t index; //!< index or pixel currently referred to }; //! returns iterator to first pixel if this material iterator begin() { return std::move(iterator(*this)); } //! returns iterator past the last pixel in this material iterator end() { return std::move(iterator(*this, false)); } protected: MaterialMuSpectre & material; //!< reference to the proxied material const StrainField_t & strain_field; //!< cell's global strain field //! references to the global stress field and perhaps tangent stiffness StressFieldTup stress_tup; //! references to the internal variables InternalVariables & internals; private: }; /* ---------------------------------------------------------------------- */ template template bool MaterialMuSpectre::iterable_proxy::iterator:: operator!=(const iterator & other) const { return (this->index != other.index); } /* ---------------------------------------------------------------------- */ template template typename MaterialMuSpectre::template iterable_proxy::iterator & MaterialMuSpectre::iterable_proxy::iterator:: operator++() { this->index++; return *this; } /* ---------------------------------------------------------------------- */ template template typename MaterialMuSpectre::template iterable_proxy< NeedTgT>::iterator::value_type MaterialMuSpectre::iterable_proxy::iterator::operator*() { const Ccoord_t pixel{ this->it.material.internal_fields.get_ccoord(this->index)}; auto && strain = std::make_tuple(this->strain_map[pixel]); auto && stresses = apply( [&pixel](auto &&... stress_tgt) { return std::make_tuple(stress_tgt[pixel]...); }, this->stress_map); auto && internal = this->it.material.get_internals(); const auto id{this->index}; auto && internals = apply( [id](auto &&... internals_) { return InternalReferences{internals_[id]...}; }, internal); return std::make_tuple(std::move(strain), std::move(stresses), std::move(internals)); } } // namespace muSpectre #endif // SRC_MATERIALS_MATERIAL_MUSPECTRE_BASE_HH_ diff --git a/tests/python_binding_tests.py b/tests/python_binding_tests.py index dd7b4ef..14e4294 100755 --- a/tests/python_binding_tests.py +++ b/tests/python_binding_tests.py @@ -1,209 +1,210 @@ #!/usr/bin/env python3 """ file python_binding_tests.py @author Till Junge @date 09 Jan 2018 @brief Unit tests for python bindings @section LICENCE Copyright © 2018 Till Junge µSpectre is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3, or (at your option) any later version. µSpectre 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 Lesser General Public License along with µSpectre; see the file COPYING. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. Additional permission under GNU GPL version 3 section 7 If you modify this Program, or any covered work, by linking or combining it with proprietary FFT implementations or numerical libraries, containing parts covered by the terms of those libraries' licenses, the licensors of this Program grant you additional permission to convey the resulting work. """ import unittest import numpy as np from python_test_imports import µ from python_fft_tests import FFT_Check from python_projection_tests import * from python_material_linear_elastic3_test import MaterialLinearElastic3_Check from python_material_linear_elastic4_test import MaterialLinearElastic4_Check from python_material_linear_elastic_generic1_test import MaterialLinearElasticGeneric1_Check from python_material_linear_elastic_generic2_test import MaterialLinearElasticGeneric2_Check from python_field_tests import FieldCollection_Check from python_exact_reference_elastic_test import LinearElastic_Check from python_field_tests import FieldCollection_Check +from python_material_evaluator_test import MaterialEvaluator_Check class CellCheck(unittest.TestCase): def test_Construction(self): """ Simple check for cell constructors """ resolution = [5,7] lengths = [5.2, 8.3] formulation = µ.Formulation.small_strain try: sys = µ.Cell(resolution, lengths, formulation) mat = µ.material.MaterialLinearElastic1_2d.make(sys, "material", 210e9, .33) except Exception as err: print(err) raise err class MaterialLinearElastic1_2dCheck(unittest.TestCase): def setUp(self): self.resolution = [5,7] self.lengths = [5.2, 8.3] self.formulation = µ.Formulation.small_strain self.sys = µ.Cell(self.resolution, self.lengths, self.formulation) self.mat = µ.material.MaterialLinearElastic1_2d.make( self.sys, "material", 210e9, .33) def test_add_material(self): self.mat.add_pixel([2,1]) class SolverCheck(unittest.TestCase): def setUp(self): self.resolution = [3, 3]#[5,7] self.lengths = [3., 3.]#[5.2, 8.3] self.formulation = µ.Formulation.finite_strain self.sys = µ.Cell(self.resolution, self.lengths, self.formulation) self.hard = µ.material.MaterialLinearElastic1_2d.make( self.sys, "hard", 210e9, .33) self.soft = µ.material.MaterialLinearElastic1_2d.make( self.sys, "soft", 70e9, .33) def test_solve(self): for i, pixel in enumerate(self.sys): if i < 3: self.hard.add_pixel(pixel) else: self.soft.add_pixel(pixel) self.sys.initialise() tol = 1e-6 Del0 = np.array([[0, .1], [0, 0]]) maxiter = 100 verbose = 0 solver=µ.solvers.SolverCG(self.sys, tol, maxiter, verbose) r = µ.solvers.de_geus(self.sys, Del0, solver,tol, verbose) #print(r) class EigenStrainCheck(unittest.TestCase): def setUp(self): self.resolution = [3, 3]#[5,7] self.lengths = [3., 3.]#[5.2, 8.3] self.formulation = µ.Formulation.small_strain self.cell1 = µ.Cell(self.resolution, self.lengths, self.formulation) self.cell2 = µ.Cell(self.resolution, self.lengths, self.formulation) self.mat1 = µ.material.MaterialLinearElastic1_2d.make( self.cell1, "simple", 210e9, .33) self.mat2 = µ.material.MaterialLinearElastic2_2d.make( self.cell2, "eigen", 210e9, .33) self.mat3 = µ.material.MaterialLinearElastic2_2d.make( self.cell2, "eigen2", 120e9, .33) def test_globalisation(self): for pixel in self.cell2: self.mat2.add_pixel(pixel, np.random.rand(2,2)) loc_eigenstrain = self.mat2.collection.get_real_field("Eigenstrain").array glo_eigenstrain = self.cell2.get_globalised_internal_real_array("Eigenstrain") error = np.linalg.norm(loc_eigenstrain-glo_eigenstrain) self.assertEqual(error, 0) def test_globalisation_constant(self): for i, pixel in enumerate(self.cell2): if i%2 == 0: self.mat2.add_pixel(pixel, np.ones((2,2))) else: self.mat3.add_pixel(pixel, np.ones((2,2))) glo_eigenstrain = self.cell2.get_globalised_internal_real_array("Eigenstrain") error = np.linalg.norm(glo_eigenstrain-1) self.assertEqual(error, 0) def test_solve(self): verbose_test = False if verbose_test: print("start test_solve") grad = np.array([[1.1, .2], [ .3, 1.5]]) gl_strain = -0.5*(grad.T.dot(grad) - np.eye(2)) gl_strain = -0.5*(grad.T + grad - 2*np.eye(2)) grad = -gl_strain if verbose_test: print("grad =\n{}\ngl_strain =\n{}".format(grad, gl_strain)) for i, pixel in enumerate(self.cell1): self.mat1.add_pixel(pixel) self.mat2.add_pixel(pixel, gl_strain) self.cell1.initialise() self.cell2.initialise() tol = 1e-6 Del0_1 = grad Del0_2 = np.zeros_like(grad) maxiter = 2 verbose = 0 def solve(cell, grad): solver=µ.solvers.SolverCG(cell, tol, maxiter, verbose) r = µ.solvers.newton_cg(cell, grad, solver, tol, tol, verbose) return r results = [solve(cell, del0) for (cell, del0) in zip((self.cell1, self.cell2), (Del0_1, Del0_2))] P1 = results[0].stress P2 = results[1].stress error = np.linalg.norm(P1-P2)/np.linalg.norm(.5*(P1+P2)) if verbose_test: print("cell 1, no eigenstrain") print("P1:\n{}".format(P1[:,0])) print("F1:\n{}".format(results[0].grad[:,0])) print("cell 2, with eigenstrain") print("P2:\n{}".format(P2[:,0])) print("F2:\n{}".format(results[1].grad[:,0])) print("end test_solve") self.assertLess(error, tol) if __name__ == '__main__': unittest.main() diff --git a/tests/python_material_evaluator_test.py b/tests/python_material_evaluator_test.py new file mode 100644 index 0000000..e0fef76 --- /dev/null +++ b/tests/python_material_evaluator_test.py @@ -0,0 +1,94 @@ +#!/usr/bin/env python3 +# -*- coding:utf-8 -*- +""" +@file python_material_evaluator_test.py + +@author Till Junge + +@date 17 Jan 2019 + +@brief tests the python bindings of the material evaluator + +Copyright © 2019 Till Junge + +µSpectre is free software; you can redistribute it and/or +modify it under the terms of the GNU General Lesser Public License as +published by the Free Software Foundation, either version 3, or (at +your option) any later version. + +µSpectre 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 Lesser General Public License +along with µSpectre; see the file COPYING. If not, write to the +Free Software Foundation, Inc., 59 Temple Place - Suite 330, +Boston, MA 02111-1307, USA. + +Additional permission under GNU GPL version 3 section 7 + +If you modify this Program, or any covered work, by linking or combining it +with proprietary FFT implementations or numerical libraries, containing parts +covered by the terms of those libraries' licenses, the licensors of this +Program grant you additional permission to convey the resulting work. +""" + + +import unittest +import numpy as np + +from python_test_imports import µ +LinMat = µ.material.MaterialLinearElastic1_2d + +def errfun(a, b): + return np.linalg.norm(a-b)/np.linalg.norm(a+b) + +class MaterialEvaluator_Check(unittest.TestCase): + def test_linear_elasticity(self): + young, poisson = 210e9, .33 + material, evaluator = LinMat.make_evaluator(young, poisson) + material.add_pixel([0,0]) + stress = evaluator.evaluate_stress(np.array([[1., 0],[0, 1]]), + µ.Formulation.finite_strain) + stress=stress.copy() + self.assertEqual(np.linalg.norm(stress), 0) + + stress2, tangent = evaluator.evaluate_stress_tangent( + np.array([[1., .0],[0, 1.0]]), + µ.Formulation.finite_strain) + + tangent = tangent.copy() + + self.assertEqual(np.linalg.norm(stress-stress2), 0) + + num_tangent = evaluator.estimate_tangent( + np.array([[1., .0],[0, 1.0]]), + µ.Formulation.finite_strain, + 1e-6, µ.FiniteDiff.centred) + tol = 1e-8 + err = errfun(num_tangent, tangent) + if not err < tol: + print("tangent:\n{}".format(tangent)) + print("num_tangent:\n{}".format(num_tangent)) + self.assertLess(err, tol) + num_tangent = evaluator.estimate_tangent( + np.array([[1., .0],[0, 1.0]]), + µ.Formulation.finite_strain, + 1e-6) + + numlin_tangent= evaluator.estimate_tangent( + np.array([[1, .0],[0, 1.0]]), + µ.Formulation.small_strain, + 1e-6, µ.FiniteDiff.centred) + + lin_stress, lin_tangent = evaluator.evaluate_stress_tangent( + np.array([[0, .0],[0, .0]]), + µ.Formulation.small_strain) + + + err = errfun(numlin_tangent, lin_tangent) + self.assertLess(err, tol) + +if __name__ == '__main__': + unittest.main() diff --git a/tests/python_test_imports.py b/tests/python_test_imports.py index 57f4351..be9192b 100644 --- a/tests/python_test_imports.py +++ b/tests/python_test_imports.py @@ -1,52 +1,53 @@ #!/usr/bin/env python3 # -*- coding:utf-8 -*- """ file python_test_imports.py @author Till Junge @date 18 Jan 2018 @brief prepares sys.path to load muSpectre @section LICENSE Copyright © 2018 Till Junge µSpectre is free software; you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3, or (at your option) any later version. µSpectre 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 Lesser General Public License along with µSpectre; see the file COPYING. If not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. Additional permission under GNU GPL version 3 section 7 If you modify this Program, or any covered work, by linking or combining it with proprietary FFT implementations or numerical libraries, containing parts covered by the terms of those libraries' licenses, the licensors of this Program grant you additional permission to convey the resulting work. """ import sys import os # Default path of the library sys.path.insert(0, os.path.join(os.getcwd(), "language_bindings/python")) # Path of the library when compiling with Xcode sys.path.insert(0, os.path.join(os.getcwd(), "language_bindings/python/Debug")) try: import muSpectre as µ + import muSpectre except ImportError as err: print(err) sys.exit(-1) diff --git a/tests/test_material_evaluator.cc b/tests/test_material_evaluator.cc new file mode 100644 index 0000000..44630ea --- /dev/null +++ b/tests/test_material_evaluator.cc @@ -0,0 +1,288 @@ +/** + * @file test_material_evaluator.cc + * + * @author Till Junge + * + * @date 13 Jan 2019 + * + * @brief tests for the material evaluator mechanism + * + * Copyright © 2019 Till Junge + * + * µSpectre is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License as + * published by the Free Software Foundation, either version 3, or (at + * your option) any later version. + * + * µSpectre 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 Lesser General Public License + * along with µSpectre; see the file COPYING. If not, write to the + * Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + * + * Additional permission under GNU GPL version 3 section 7 + * + * If you modify this Program, or any covered work, by linking or combining it + * with proprietary FFT implementations or numerical libraries, containing parts + * covered by the terms of those libraries' licenses, the licensors of this + * Program grant you additional permission to convey the resulting work. + */ + +#include "tests.hh" +#include "common/T4_map_proxy.hh" +#include "materials/material_linear_elastic2.hh" +#include "materials/material_evaluator.hh" + +#include "Eigen/Dense" + +namespace muSpectre { + + BOOST_AUTO_TEST_SUITE(material_evaluator_tests); + + /* ---------------------------------------------------------------------- */ + BOOST_AUTO_TEST_CASE(without_per_pixel_data) { + using Mat_t = MaterialLinearElastic1; + + constexpr Real Young{210e9}; + constexpr Real Poisson{.33}; + + auto mat_eval = Mat_t::make_evaluator(Young, Poisson); + auto & mat = *std::get<0>(mat_eval); + auto & evaluator = std::get<1>(mat_eval); + + using T2_t = Eigen::Matrix; + using T4_t = T4Mat; + const T2_t F{(T2_t::Random() - (T2_t::Ones() * .5)) * 1e-4 + + T2_t::Identity()}; + const T2_t eps{ + .5 * ((F - T2_t::Identity()) + (F - T2_t::Identity()).transpose())}; + + /* + * at this point, the evaluator has been created, but the underlying + * material still has zero pixels. Evaluation is not yet possible, and + * trying to do so has to fail with an explicit error message + */ + BOOST_CHECK_THROW(evaluator.evaluate_stress(eps, Formulation::small_strain), + std::runtime_error); + + mat.add_pixel({}); + + const T2_t sigma{evaluator.evaluate_stress(eps, Formulation::small_strain)}; + const T2_t P{evaluator.evaluate_stress(F, Formulation::finite_strain)}; + + auto J{F.determinant()}; + auto P_reconstruct{J * sigma * F.inverse().transpose()}; + + auto error_comp{[](const auto & a, const auto & b) { + return (a - b).norm() / (a + b).norm(); + }}; + auto error{error_comp(P, P_reconstruct)}; + + constexpr Real small_strain_tol{1e-3}; + if (not(error <= small_strain_tol)) { + std::cout << "F =" << std::endl << F << std::endl; + std::cout << "ε =" << std::endl << eps << std::endl; + std::cout << "P =" << std::endl << P << std::endl; + std::cout << "σ =" << std::endl << sigma << std::endl; + std::cout << "P_reconstructed =" << std::endl + << P_reconstruct << std::endl; + } + + BOOST_CHECK_LE(error, small_strain_tol); + + T2_t sigma2, P2; + T4_t C, K; + + std::tie(sigma2, C) = + evaluator.evaluate_stress_tangent(eps, Formulation::small_strain); + std::tie(P2, K) = + evaluator.evaluate_stress_tangent(F, Formulation::finite_strain); + + error = error_comp(sigma2, sigma); + BOOST_CHECK_LE(error, tol); + error = error_comp(P2, P); + BOOST_CHECK_LE(error, tol); + + error = error_comp(C, K); + if (not(error <= small_strain_tol)) { + std::cout << "F =" << std::endl << F << std::endl; + std::cout << "ε =" << std::endl << eps << std::endl; + std::cout << "P =" << std::endl << P << std::endl; + std::cout << "σ =" << std::endl << sigma << std::endl; + std::cout << "K =" << std::endl << K << std::endl; + std::cout << "C =" << std::endl << C << std::endl; + } + BOOST_CHECK_LE(error, small_strain_tol); + + + mat.add_pixel({1}); + /* + * Now, the material has two pixels, and evaluating it would be ambiguous. + * It should fail with an explicit error message + */ + BOOST_CHECK_THROW(evaluator.evaluate_stress(eps, Formulation::small_strain), + std::runtime_error); + } + + /* ---------------------------------------------------------------------- */ + BOOST_AUTO_TEST_CASE(with_per_pixel_data) { + using Mat_t = MaterialLinearElastic2; + + constexpr Real Young{210e9}; + constexpr Real Poisson{.33}; + + auto mat_eval{Mat_t::make_evaluator(Young, Poisson)}; + auto & mat{*std::get<0>(mat_eval)}; + auto & evaluator{std::get<1>(mat_eval)}; + + + using T2_t = Eigen::Matrix; + using T4_t = T4Mat; + const T2_t F{(T2_t::Random() - (T2_t::Ones() * .5)) * 1e-4 + + T2_t::Identity()}; + const T2_t eps{ + .5 * ((F - T2_t::Identity()) + (F - T2_t::Identity()).transpose())}; + + BOOST_CHECK_THROW(evaluator.evaluate_stress(eps, Formulation::small_strain), + std::runtime_error); + + T2_t eigen_strain{[](auto x) { + return 1e-4 * (x + x.transpose()); + }(T2_t::Random() - T2_t::Ones() * .5)}; + + mat.add_pixel({}, eigen_strain); + + const T2_t sigma{evaluator.evaluate_stress(eps, Formulation::small_strain)}; + const T2_t P{evaluator.evaluate_stress(F, Formulation::finite_strain)}; + + auto J{F.determinant()}; + auto P_reconstruct{J * sigma * F.inverse().transpose()}; + + auto error_comp{[](const auto & a, const auto & b) { + return (a - b).norm() / (a + b).norm(); + }}; + auto error{error_comp(P, P_reconstruct)}; + + constexpr Real small_strain_tol{1e-3}; + if (not(error <= small_strain_tol)) { + std::cout << "F =" << std::endl << F << std::endl; + std::cout << "ε =" << std::endl << eps << std::endl; + std::cout << "P =" << std::endl << P << std::endl; + std::cout << "σ =" << std::endl << sigma << std::endl; + std::cout << "P_reconstructed =" << std::endl + << P_reconstruct << std::endl; + } + + BOOST_CHECK_LE(error, small_strain_tol); + + T2_t sigma2, P2; + T4_t C, K; + + std::tie(sigma2, C) = + evaluator.evaluate_stress_tangent(eps, Formulation::small_strain); + std::tie(P2, K) = + evaluator.evaluate_stress_tangent(F, Formulation::finite_strain); + + error = error_comp(sigma2, sigma); + BOOST_CHECK_LE(error, tol); + error = error_comp(P2, P); + BOOST_CHECK_LE(error, tol); + + error = error_comp(C, K); + if (not(error <= small_strain_tol)) { + std::cout << "F =" << std::endl << F << std::endl; + std::cout << "ε =" << std::endl << eps << std::endl; + std::cout << "P =" << std::endl << P << std::endl; + std::cout << "σ =" << std::endl << sigma << std::endl; + std::cout << "K =" << std::endl << K << std::endl; + std::cout << "C =" << std::endl << C << std::endl; + } + BOOST_CHECK_LE(error, small_strain_tol); + } + + /* ---------------------------------------------------------------------- */ + BOOST_AUTO_TEST_CASE(tangent_estimation) { + using Mat_t = MaterialLinearElastic1; + + constexpr Real Young{210e9}; + constexpr Real Poisson{.33}; + + auto mat_eval = Mat_t::make_evaluator(Young, Poisson); + auto & mat = *std::get<0>(mat_eval); + auto & evaluator = std::get<1>(mat_eval); + + using T2_t = Eigen::Matrix; + using T4_t = T4Mat; + const T2_t F{(T2_t::Random() - (T2_t::Ones() * .5)) * 1e-4 + + T2_t::Identity()}; + const T2_t eps{ + .5 * ((F - T2_t::Identity()) + (F - T2_t::Identity()).transpose())}; + + BOOST_CHECK_THROW(evaluator.evaluate_stress(eps, Formulation::small_strain), + std::runtime_error); + + mat.add_pixel({}); + + T2_t sigma, P; + T4_t C, K; + + std::tie(sigma, C) = + evaluator.evaluate_stress_tangent(eps, Formulation::small_strain); + std::tie(P, K) = + evaluator.evaluate_stress_tangent(F, Formulation::finite_strain); + + constexpr Real linear_step{1.}; + constexpr Real nonlin_step{1.e-6}; + T4_t C_estim{evaluator.estimate_tangent(eps, Formulation::small_strain, + linear_step)}; + T4_t K_estim{evaluator.estimate_tangent(F, Formulation::finite_strain, + nonlin_step)}; + + auto error_comp{[](const auto & a, const auto & b) { + return (a - b).norm() / (a + b).norm(); + }}; + + constexpr Real finite_diff_tol{1e-9}; + Real error {error_comp(K, K_estim)}; + if (not(error <= finite_diff_tol)) { + std::cout << "K =" << std::endl << K << std::endl; + std::cout << "K_estim =" << std::endl << K_estim << std::endl; + } + BOOST_CHECK_LE(error, finite_diff_tol); + + error = error_comp(C, C_estim); + if (not(error <= tol)) { + std::cout << "centred difference:" << std::endl; + std::cout << "C =" << std::endl << C << std::endl; + std::cout << "C_estim =" << std::endl << C_estim << std::endl; + } + BOOST_CHECK_LE(error, tol); + + C_estim = evaluator.estimate_tangent(eps, Formulation::small_strain, + linear_step, FiniteDiff::forward); + error = error_comp(C, C_estim); + if (not(error <= tol)) { + std::cout << "forward difference:" << std::endl; + std::cout << "C =" << std::endl << C << std::endl; + std::cout << "C_estim =" << std::endl << C_estim << std::endl; + } + BOOST_CHECK_LE(error, tol); + + C_estim = evaluator.estimate_tangent(eps, Formulation::small_strain, + linear_step, FiniteDiff::backward); + error = error_comp(C, C_estim); + if (not(error <= tol)) { + std::cout << "backward difference:" << std::endl; + std::cout << "C =" << std::endl << C << std::endl; + std::cout << "C_estim =" << std::endl << C_estim << std::endl; + } + BOOST_CHECK_LE(error, tol); +} + + BOOST_AUTO_TEST_SUITE_END(); + +} // namespace muSpectre