diff --git a/language_bindings/python/bind_py_material.cc b/language_bindings/python/bind_py_material.cc
index 667feda..0e09ef5 100644
--- a/language_bindings/python/bind_py_material.cc
+++ b/language_bindings/python/bind_py_material.cc
@@ -1,190 +1,252 @@
/**
* @file bind_py_material.cc
*
* @author Till Junge <till.junge@epfl.ch>
*
* @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 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 GNU Emacs; 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 "materials/material_linear_elastic2.hh"
#include "materials/material_linear_elastic3.hh"
#include "materials/material_linear_elastic4.hh"
+#include "materials/material_base.hh"
#include "cell/cell_base.hh"
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/eigen.h>
#include <sstream>
#include <string>
using namespace muSpectre;
namespace py = pybind11;
using namespace pybind11::literals;
+
+
+template <Dim_t dim>
+void add_material_base_helper(py::module & mod) {
+ std::stringstream name_stream{};
+ name_stream << "MaterialBase" << dim << 'd';
+ const auto name {name_stream.str()};
+ using Mat_t = MaterialBase<dim, dim>;
+ py::class_<Mat_t>(mod, name.c_str());
+}
+
/**
* python binding for the optionally objective form of Hooke's law
*/
template <Dim_t dim>
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 = MaterialLinearElastic1<dim, dim>;
using Sys_t = CellBase<dim, dim>;
- py::class_<Mat_t>(mod, name.c_str())
+ using MatBase_t = MaterialBase<dim,dim>;
+ py::class_<Mat_t, MatBase_t>(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, Ccoord_t<dim> pix) {
mat.add_pixel(pix);},
"pixel"_a)
.def("add_pixel_split",
[] (Mat_t & mat, Ccoord_t<dim> pix, Real ratio) {
mat.add_pixel_split(pix, ratio);},
"pixel"_a,
"ratio"_a)
.def("size", &Mat_t::size);
}
template <Dim_t dim>
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 = MaterialLinearElastic2<dim, dim>;
using Sys_t = CellBase<dim, dim>;
- py::class_<Mat_t>(mod, name.c_str())
+ using MatBase_t = MaterialBase<dim,dim>;
+ py::class_<Mat_t, MatBase_t>(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, Ccoord_t<dim> pix, py::EigenDRef<Eigen::ArrayXXd>& eig) {
Eigen::Matrix<Real, dim, dim> eig_strain{eig};
mat.add_pixel(pix, eig_strain);},
"pixel"_a,
"eigenstrain"_a)
.def("add_pixel_split",
[] (Mat_t & mat, Ccoord_t<dim> pix, Real ratio, py::EigenDRef<Eigen::ArrayXXd>& eig) {
Eigen::Matrix<Real, dim, dim> eig_strain{eig};
mat.add_pixel_split(pix, ratio, eig_strain);},
"pixel"_a,
"ratio"_a,
"eigenstrain"_a)
.def("size", &Mat_t::size);
}
template <Dim_t dim>
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 = MaterialLinearElastic3<dim, dim>;
using Sys_t = CellBase<dim, dim>;
- py::class_<Mat_t>(mod, name.c_str())
+ using MatBase_t = MaterialBase<dim,dim>;
+ py::class_<Mat_t, MatBase_t>(mod, name.c_str())
.def(py::init<std::string>(), "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, Ccoord_t<dim> pix, Real Young, Real Poisson) {
mat.add_pixel(pix, Young, Poisson);},
"pixel"_a,
"Young"_a,
"Poisson"_a)
.def("add_pixel_split",
[] (Mat_t & mat, Ccoord_t<dim> pix, Real ratio, Real Young, Real Poisson) {
mat.add_pixel_split(pix, ratio, Young, Poisson);},
"pixel"_a,
"ratio"_a,
"Young"_a,
"Poisson"_a)
.def("size", &Mat_t::size);
}
template <Dim_t dim>
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 = MaterialLinearElastic4<dim, dim>;
using Sys_t = CellBase<dim, dim>;
py::class_<Mat_t>(mod, name.c_str())
.def(py::init<std::string>(), "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, Ccoord_t<dim> pix, Real Young, Real Poisson) {
mat.add_pixel(pix, Young, Poisson);},
"pixel"_a,
"Young"_a,
"Poisson"_a)
.def("add_pixel_split",
[] (Mat_t & mat, Ccoord_t<dim> pix, Real ratio, Real Young, Real Poisson) {
mat.add_pixel_split(pix, ratio, Young, Poisson);},
"pixel"_a,
"ratio"_a,
"Young"_a,
"Poisson"_a)
.def("size", &Mat_t::size);
}
-
+template <Dim_t Dim>
+class PyMaterialBase : public MaterialBase<Dim, Dim> {
+public:
+ /* Inherit the constructors */
+ using Parent = MaterialBase<Dim, Dim>;
+ 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) */
+ );
+ }
+ virtual void compute_stresses(const typename Parent::StrainField_t & F,
+ typename Parent::StressField_t & P,
+ 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
+ );
+ }
+
+ virtual void compute_stresses_tangent(const typename Parent::StrainField_t & F,
+ typename Parent::StressField_t & P,
+ typename Parent::TangentField_t & K,
+ Formulation form) override {
+ PYBIND11_OVERLOAD_PURE
+ (void, /* Return type */
+ Parent, /* Parent class */
+ compute_stresses_tangent, /* Name of function in C++ (must match Python name) */
+ F,P,K, form
+ );
+ }
+
+};
template <Dim_t dim>
void add_material_helper(py::module & mod) {
+ add_material_base_helper<dim>(mod);
add_material_linear_elastic1_helper<dim>(mod);
add_material_linear_elastic2_helper<dim>(mod);
add_material_linear_elastic3_helper<dim>(mod);
add_material_linear_elastic4_helper<dim>(mod);
}
void add_material(py::module & mod) {
auto material{mod.def_submodule("material")};
material.doc() = "bindings for constitutive laws";
add_material_helper<twoD >(material);
add_material_helper<threeD>(material);
}
diff --git a/src/common/intersection_octree.hh b/src/common/intersection_octree.hh
index b9fac69..4465a54 100644
--- a/src/common/intersection_octree.hh
+++ b/src/common/intersection_octree.hh
@@ -1,145 +1,145 @@
/**
* @file intersection_octree.hh
*
* @author Ali Falsafi <ali.falsafi@epfl.ch>
*
* @date May 2018
*
* @brief common operations on pixel addressing
*
* Copyright © 2018 Ali Falsafi
*
* µ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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#ifndef INTERSECTION_OCTREE_H
#define INTERSECTION_OCTREE_H
#include <vector>
#include "common/ccoord_operations.hh"
#include "common/common.hh"
#include "cell/cell_base.hh"
#include "materials/material_base.hh"
#include "common/intersection_volume_calculator.hh"
namespace muSpectre {
template<Dim_t DimS>
class RootNode;
template<Dim_t DimS>
class Node
{
public:
using Rcoord = Rcoord_t<DimS>; //!< physical coordinates type
using Ccoord = Ccoord_t<DimS>; //!< cell coordinates type
using RootNode_t = RootNode<DimS>;
//! Default constructor
Node() = delete;
// Construct by origin, lenghts, and depth
Node(const Rcoord &new_origin, const Ccoord &new_lenghts,
int depth, RootNode_t& root, bool is_root );
// Construct by cell
// Node(const CellBase<DimS, DimS>& cell, RootNode_t& root);
//! Copy constructor
Node(const Node &other) = delete;
//! Move constructor
Node(Node &&other) = default;
//! Destructor
- ~Node() = default;
+ virtual ~Node() = default;
//This function checks the status of the node and orders its devision into
//smaller nodes or asssign material to it
virtual void check_node();
//This function gives the ratio of the node which happens to be inside the precipitate
//and assign materila to it if node is a pixel or divide it furhter if it is not
void split_node(Real ratio);
//this function constructs children of a node
void divide_node();
protected:
////////
RootNode_t& root_node;
Rcoord origin, Rlengths{};
Ccoord Clengths{};
int depth;
bool is_pixel;
int children_no;
std::vector<Node> children{};
};
template<Dim_t DimS>
class RootNode: public Node<DimS>
{
friend class Node<DimS>;
public:
using Rcoord = Rcoord_t<DimS>; //!< physical coordinates type
using Ccoord = Ccoord_t<DimS>; //!< cell coordinates type
using Parent = Node<DimS>; //!< base class
//!Default Constructor
RootNode() = delete ;
//! Constructing a root node for a cell and a preticipate inside that cell
RootNode(CellBase<DimS, DimS>& cell,
std::vector<Rcoord> vert_precipitate);
//! Copy constructor
RootNode(const RootNode &other) = delete;
//! Move constructor
RootNode(RootNode &&other) = default;
//! Destructor
~RootNode() = default;
//returns the pixels which have intersection raio with the preipitate
std::vector<Ccoord> get_intersected_pixels ();
//return the intersection ratios of corresponding to the pixels returned by get_intersected_pixels()
std::vector<Real> get_intersection_ratios ();
//checking rootnode:
void check_root_node();
protected:
CellBase<DimS, DimS>& cell;
Rcoord cell_length, pixel_lengths;
Ccoord cell_resolution;
int max_resolution ;
int max_depth;
std::vector<Rcoord> precipitate_vertices{};
std::vector<Ccoord> intersected_pixels{};
std::vector<Real> intersection_ratios{};
};
} //muSpectre
#endif /* INTERSECTION_OCTREE_H */
diff --git a/src/materials/material_muSpectre_base.hh b/src/materials/material_muSpectre_base.hh
index 932f2b8..2eb3036 100644
--- a/src/materials/material_muSpectre_base.hh
+++ b/src/materials/material_muSpectre_base.hh
@@ -1,911 +1,911 @@
/**
* @file material_muSpectre_base.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @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 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 GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#ifndef MATERIAL_MUSPECTRE_BASE_H
#define MATERIAL_MUSPECTRE_BASE_H
#include "common/common.hh"
#include "materials/material_base.hh"
#include "materials/materials_toolbox.hh"
#include "common/field_collection.hh"
#include "common/field.hh"
#include "common//utilities.hh"
#include <tuple>
#include <type_traits>
#include <iterator>
#include <stdexcept>
namespace muSpectre {
// Forward declaration for factory function
template <Dim_t DimS, Dim_t DimM>
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 <class Material>
struct MaterialMuSpectre_traits {
};
template <class Material, Dim_t DimS, Dim_t DimM>
class MaterialMuSpectre;
/**
* Base class for most convenient implementation of materials
*/
template <class Material, Dim_t DimS, Dim_t DimM>
class MaterialMuSpectre: public MaterialBase<DimS, DimM>
{
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<DimS, DimM>; //!< 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;
using Ccoord = Ccoord_t<DimS>; //!< cell coordinates type
//! traits for the CRTP subclass
using traits = MaterialMuSpectre_traits<Material>;
//! Default constructor
MaterialMuSpectre() = delete;
//! Construct by name
MaterialMuSpectre(std::string name);
//! Copy constructor
MaterialMuSpectre(const MaterialMuSpectre &other) = delete;
//! Move constructor
MaterialMuSpectre(MaterialMuSpectre &&other) = delete;
//! Destructor
virtual ~MaterialMuSpectre() = default;
//! Factory
template <class... ConstructorArgs>
static Material & make(CellBase<DimS, DimM> & cell,
ConstructorArgs &&... args);
//! Copy assignment operator
MaterialMuSpectre& operator=(const MaterialMuSpectre &other) = delete;
//! Move assignment operator
MaterialMuSpectre& operator=(MaterialMuSpectre &&other) = delete;
//! allocate memory, etc
virtual void initialise() override;
// this fucntion is speicalized to assign partial material to a pixel
//void add_pixel_split(const Ccoord & local_ccoord, Real ratio = 1.0);
template<class ...InternalArgs>
inline void add_pixel_split(const Ccoord_t<DimS> & pixel,
Real ratio,
InternalArgs ... args);
inline void add_pixel_split(const Ccoord_t<DimS> & pixel,
- Real ratio);
+ Real ratio) override;
void add_split_pixels_precipitate(std::vector<Ccoord_t<DimS>> intersected_precipitate,
std::vector<Real> intersection_ratios);
//! computes stress
using Parent::compute_stresses;
using Parent::compute_stresses_tangent;
virtual void compute_stresses(const StrainField_t & F,
StressField_t & P,
Formulation form,
SplittedCell is_cell_splitted) override;
//! computes stress and tangent modulus
virtual void compute_stresses_tangent(const StrainField_t & F,
StressField_t & P,
TangentField_t & K,
Formulation form,
SplittedCell is_cell_splitted)override;
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 <Formulation Form, SplittedCell is_cell_splitted>
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 <Formulation Form, SplittedCell is_cell_splitted>
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<NeedTangent need_tgt = NeedTangent::no>
class iterable_proxy;
/**
* inheriting classes with internal variables need to overload this function
*/
typename traits::InternalVariables& get_internals() {
return static_cast<Material&>(*this).get_internals();}
bool is_initialised{false}; //!< to handle double initialisation right
private:
};
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
MaterialMuSpectre<Material, DimS, DimM>::
MaterialMuSpectre(std::string name)
:Parent(name){
using stress_compatible = typename traits::StressMap_t::
template is_compatible<StressField_t>;
using strain_compatible = typename traits::StrainMap_t::
template is_compatible<StrainField_t>;
using tangent_compatible = typename traits::TangentMap_t::
template is_compatible<TangentField_t>;
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 <class Material, Dim_t DimS, Dim_t DimM>
template <class... ConstructorArgs>
Material & MaterialMuSpectre<Material, DimS, DimM>::
make(CellBase<DimS, DimM> & cell,
ConstructorArgs && ... args) {
auto mat = std::make_unique<Material>(args...);
auto & mat_ref = *mat;
cell.add_material(std::move(mat));
return mat_ref;
}
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
void MaterialMuSpectre<Material, DimS, DimM>::
add_pixel_split(const Ccoord &local_ccoord, Real ratio) {
auto & this_mat = static_cast<Material&>(*this);
this_mat.add_pixel(local_ccoord);
this->assigned_ratio.push_back(ratio);
}
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
template<class ...InternalArgs>
void MaterialMuSpectre<Material, DimS, DimM>::
add_pixel_split(const Ccoord_t<DimS> & pixel,
Real ratio,
InternalArgs ... args){
auto & this_mat = static_cast<Material&>(*this);
this_mat.add_pixel(pixel, args...);
this->assigned_ratio.push_back(ratio);
}
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
void MaterialMuSpectre<Material, DimS, DimM>::
add_split_pixels_precipitate(std::vector<Ccoord_t<DimS>> intersected_precipitate,
std::vector<Real> intersection_ratios){
//assign precipitate materials:
for (auto && tup: akantu::zip(intersected_precipitate, intersection_ratios) ){
auto pix = std::get<0> (tup);
auto ratio = std::get<1> (tup);
this->add_pixel_split(pix,ratio);
}
}
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
void MaterialMuSpectre<Material, DimS, DimM>::
initialise() {
if (!this->is_initialised) {
this->internal_fields.initialise();
this->is_initialised = true;
}
}
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
void MaterialMuSpectre<Material, DimS, DimM>::
compute_stresses(const StrainField_t &F, StressField_t &P,
Formulation form, SplittedCell is_cell_splitted) {
switch (form) {
case Formulation::finite_strain: {
switch (is_cell_splitted){
case (SplittedCell::no): {
this->compute_stresses_worker<Formulation::finite_strain, SplittedCell::no>(F, P);
break;
}
case (SplittedCell::yes): {
this->compute_stresses_worker<Formulation::finite_strain, SplittedCell::yes>(F, P);
break;
}
}
break;
}
case Formulation::small_strain: {
switch (is_cell_splitted){
case (SplittedCell::no): {
this->compute_stresses_worker<Formulation::small_strain, SplittedCell::no>(F, P);
break;
}
case (SplittedCell::yes): {
this->compute_stresses_worker<Formulation::small_strain, SplittedCell::yes>(F, P);
break;
}
}
break;
}
default:
throw std::runtime_error("Unknown formulation");
break;
}
}
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
void MaterialMuSpectre<Material, DimS, DimM>::
compute_stresses_tangent(const StrainField_t & F, StressField_t & P,
TangentField_t & K,
Formulation form, SplittedCell is_cell_splitted) {
switch (form) {
case Formulation::finite_strain: {
switch (is_cell_splitted){
case (SplittedCell::no): {
this->compute_stresses_worker<Formulation::finite_strain, SplittedCell::no>(F, P, K);
break;
}
case (SplittedCell::yes): {
this->compute_stresses_worker<Formulation::finite_strain, SplittedCell::yes>(F, P, K);
break;
}
}
break;
}
case Formulation::small_strain: {
switch (is_cell_splitted){
case (SplittedCell::no): {
this->compute_stresses_worker<Formulation::small_strain, SplittedCell::no>(F, P, K);
break;
}
case (SplittedCell::yes): {
this->compute_stresses_worker<Formulation::small_strain, SplittedCell::yes>(F, P, K);
break;
}
}
break;
}
default:
throw std::runtime_error("Unknown formulation");
break;
}
}
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
template <Formulation Form, SplittedCell is_cell_splitted>
void MaterialMuSpectre<Material, DimS, DimM>::
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<typename traits::StrainMap_t::reference>;
using Stresses_t = std::tuple<typename traits::StressMap_t::reference,
typename traits::TangentMap_t::reference>;
auto constitutive_law_small_strain = [this]
(Strains_t Strains, Stresses_t Stresses, const Real & ratio, 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<Material&>(*this);
// Transformation gradient is first in the strains tuple
auto & grad = std::get<0>(Strains);
auto && strain = MatTB::convert_strain<stored_strain_m, expected_strain_m>(grad);
// return value contains a tuple of rvalue_refs to both stress and tangent moduli
// for the case that cells do not contain any splitt cells
auto non_split_pixel = [&strain, &this_mat, ratio, &Stresses, &internal_variables](){
Stresses =
apply([&strain, &this_mat] (auto && ... internals) {
return
this_mat.evaluate_stress_tangent(std::move(strain),
internals...);},
internal_variables);
};
// for the case that cells contain splitt cells
auto split_pixel = [&strain, &this_mat, ratio, &Stresses, &internal_variables](){
auto stress_tgt_contributions =
apply([&strain, &this_mat] (auto && ... internals) {
return
this_mat.evaluate_stress_tangent(std::move(strain),
internals...);},
internal_variables);
auto && stress = std::get<0>(Stresses);
auto && tangent = std::get<1>(Stresses);
stress += ratio * std::get<0>(stress_tgt_contributions);
tangent += ratio * std::get<1>(stress_tgt_contributions);
};
switch (is_cell_splitted){
case SplittedCell::no : {
non_split_pixel();
break;
}
case SplittedCell::yes : {
split_pixel();
break;
}
}
};
auto constitutive_law_finite_strain = [this]
(Strains_t Strains, Stresses_t Stresses, const Real & ratio, 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<Material&>(*this);
// Transformation gradient is first in the strains tuple
auto & grad = std::get<0>(Strains);
auto && strain = MatTB::convert_strain<stored_strain_m, expected_strain_m>(grad);
// return value contains a tuple of rvalue_refs to both stress
// and tangent moduli
// for the case that cells do not contain splitt cells
auto non_split_pixel = [&strain, &this_mat, ratio, &Stresses, &internal_variables, &grad](){
auto stress_tgt =
apply([&strain, &this_mat] (auto && ... internals) {
return
this_mat.evaluate_stress_tangent(std::move(strain),
internals...);},
internal_variables);
Stresses = MatTB::PK1_stress<traits::stress_measure, traits::strain_measure>
(std::move(grad),
std::move(std::get<0>(stress_tgt)),
std::move( std::get<1>(stress_tgt)));
};
// for the case that cells contain splitt cells
auto split_pixel = [&strain, &this_mat, ratio, &Stresses, &internal_variables, &grad](){
auto stress_tgt =
apply([&strain, &this_mat] (auto && ... internals) {
return
this_mat.evaluate_stress_tangent(std::move(strain),
internals...);},
internal_variables);
auto && stress_tgt_contributions = MatTB::PK1_stress<traits::stress_measure, traits::strain_measure>
(std::move(grad),
std::move(std::get<0>(stress_tgt)),
std::move(std::get<1>(stress_tgt)));
auto && stress = std::get<0>(Stresses);
auto && tangent = std::get<1>(Stresses);
stress += ratio * std::get<0>(stress_tgt_contributions);
tangent += ratio * std::get<1>(stress_tgt_contributions);
};
switch (is_cell_splitted){
case SplittedCell::no : {
non_split_pixel();
break;
}
case SplittedCell::yes : {
split_pixel();
break;
}
}
};
iterable_proxy<NeedTangent::yes> 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<typename traits::StrainMap_t::reference,
std::remove_reference_t<
decltype(std::get<0>(std::get<0>(arglist)))>>::value,
"Type mismatch for strain reference, check iterator "
"value_type");
static_assert(std::is_same<typename traits::StressMap_t::reference,
std::remove_reference_t<
decltype(std::get<0>(std::get<1>(arglist)))>>::value,
"Type mismatch for stress reference, check iterator"
"value_type");
static_assert(std::is_same<typename traits::TangentMap_t::reference,
std::remove_reference_t<
decltype(std::get<1>(std::get<1>(arglist)))>>::value,
"Type mismatch for tangent reference, check iterator"
"value_type");
static_assert(std::is_same<Real,
std::remove_reference_t<
decltype(std::get<2>(arglist))>>::value,
"Type mismatch for ratio reference, expected a real number");
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 <class Material, Dim_t DimS, Dim_t DimM>
template <Formulation Form, SplittedCell is_cell_splitted>
void MaterialMuSpectre<Material, DimS, DimM>::
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<typename traits::StrainMap_t::reference>;
using Stresses_t = std::tuple<typename traits::StressMap_t::reference>;
auto constitutive_law_small_strain = [this]
(Strains_t Strains, Stresses_t Stresses,const Real & ratio, 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<Material&>(*this);
// Transformation gradient is first in the strains tuple
auto & grad = std::get<0>(Strains);
auto && strain = MatTB::convert_strain<stored_strain_m, expected_strain_m>(grad);
// return value contains a tuple of rvalue_refs to both stress and tangent moduli
auto & sigma = std::get<0>(Stresses);
//for the case that cell does not contain any splitted pixel
auto non_split_pixel = [&strain, &this_mat,ratio, &sigma, &internal_variables](){
sigma =
apply([&strain, &this_mat] (auto && ... internals) {
return
this_mat.evaluate_stress(std::move(strain),
internals...);},
internal_variables);
};
// for the case that cells contain splitt cells
auto split_pixel = [&strain, &this_mat, &ratio, &sigma, &internal_variables](){
sigma += ratio *
apply([&strain, &this_mat] (auto && ... internals) {
return
this_mat.evaluate_stress(std::move(strain),
internals...);},
internal_variables);
};
switch (is_cell_splitted){
case SplittedCell::no : {
non_split_pixel();
break;
}
case SplittedCell::yes : {
split_pixel();
break;
}
}
};
auto constitutive_law_finite_strain = [this]
(Strains_t Strains, Stresses_t && Stresses,const Real & ratio, 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<Material&>(*this);
// Transformation gradient is first in the strains tuple
auto & grad = std::get<0>(Strains);
auto && strain = MatTB::convert_strain<stored_strain_m, expected_strain_m>(grad);
// TODO: 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.
//for the case that cell does not contain any splitted pixel
auto non_split_pixel = [ &strain, &this_mat, &ratio, &Stresses, &internal_variables, &grad] (){
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<traits::stress_measure, traits::strain_measure>
(grad, stress);
};
// for the case that cells contain splitted cells
auto split_pixel = [&strain, &this_mat, ratio, &Stresses, internal_variables, &grad](){
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 += ratio * MatTB::PK1_stress<traits::stress_measure, traits::strain_measure>
(grad, stress);
};
switch (is_cell_splitted){
case SplittedCell::no : {
non_split_pixel();
break;
}
case SplittedCell::yes : {
split_pixel();
break;
}
}
};
iterable_proxy<NeedTangent::no> fields{*this, F, P};
for (auto && arglist: fields) {
/**
* arglist is a tuple of three tuples containing only Lvalue
* references (see value_type 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<typename traits::StrainMap_t::reference,
std::remove_reference_t<
decltype(std::get<0>(std::get<0>(arglist)))>>::value,
"Type mismatch for strain reference, check iterator "
"value_type");
static_assert(std::is_same<typename traits::StressMap_t::reference,
std::remove_reference_t<
decltype(std::get<0>(std::get<1>(arglist)))>>::value,
"Type mismatch for stress reference, check iterator"
"value_type");
static_assert(std::is_same<Real,
std::remove_reference_t<
decltype(std::get<2>(arglist))>>::value,
"Type mismatch for ratio reference, expected a real number");
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 <class Material, Dim_t DimS, Dim_t DimM>
template <MatTB::NeedTangent NeedTgt>
class MaterialMuSpectre<Material, DimS, DimM>::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<Material, DimS, DimM>::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<bool DoNeedTgt=(NeedTgt == NeedTangent::yes)>
iterable_proxy(MaterialMuSpectre & mat,
const StrainField_t & F,
StressField_t & P,
std::enable_if_t<DoNeedTgt, TangentField_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<bool DontNeedTgt=(NeedTgt == NeedTangent::no)>
iterable_proxy(MaterialMuSpectre & mat,
const StrainField_t & F,
std::enable_if_t<DontNeedTgt, StressField_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<StressField_t&, TangentField_t&>,
std::tuple<StressField_t&>>;
//! tuple containing a stress and possibly a tangent stiffness field map
using StressMapTup = std::conditional_t
<(NeedTgt == NeedTangent::yes),
std::tuple<StressMap_t, TangentMap_t>,
std::tuple<StressMap_t>>;
//! tuple containing a stress and possibly a tangent stiffness value ref
using Stress_tTup = std::conditional_t<(NeedTgt == NeedTangent::yes),
std::tuple<Stress_t, Tangent_t>,
std::tuple<Stress_t>>;
//! 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<InternalVariables>;
//! return type to be unpacked per pixel my the constitutive law
using value_type =
std::tuple<std::tuple<Strain_t>, Stress_tTup, Real, 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.
**/
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
private:
};
//! 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 <class Material, Dim_t DimS, Dim_t DimM>
template <MatTB::NeedTangent NeedTgt>
bool
MaterialMuSpectre<Material, DimS, DimM>::iterable_proxy<NeedTgt>::iterator::
operator!=(const iterator & other) const {
return (this->index != other.index);
}
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
template <MatTB::NeedTangent NeedTgt>
typename MaterialMuSpectre<Material, DimS, DimM>::
template iterable_proxy<NeedTgt>::
iterator &
MaterialMuSpectre<Material, DimS, DimM>::iterable_proxy<NeedTgt>::iterator::
operator++() {
this->index++;
return *this;
}
/* ---------------------------------------------------------------------- */
template <class Material, Dim_t DimS, Dim_t DimM>
template <MatTB::NeedTangent NeedTgT>
typename MaterialMuSpectre<Material, DimS, DimM>::
template iterable_proxy<NeedTgT>::iterator::
value_type
MaterialMuSpectre<Material, DimS, DimM>::iterable_proxy<NeedTgT>::iterator::
operator*() {
const Ccoord_t<DimS> pixel{
this->it.material.internal_fields.get_ccoord(this->index)};
auto && strain = std::make_tuple(this->strain_map[pixel]);
auto && ratio = this->it.material.get_assigned_ratio(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(ratio),
std::move(internals));
}
} // muSpectre
#endif /* MATERIAL_MUSPECTRE_BASE_H */