diff --git a/language_bindings/python/bind_py_material.cc b/language_bindings/python/bind_py_material.cc
index faf4d25..0997dbb 100644
--- a/language_bindings/python/bind_py_material.cc
+++ b/language_bindings/python/bind_py_material.cc
@@ -1,330 +1,332 @@
/**
* @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_anisotropic.hh"
#include "materials/material_orthotropic.hh"
#include "materials/material_base.hh"
#include "cell/cell_base.hh"
#include "common/field_collection.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 Anisotropic material
*/
template <Dim_t dim>
void add_material_anisotropic_helper(py::module & mod) {
std::stringstream name_stream{};
name_stream << "MaterialAnisotropic_" << dim << 'd';
const auto name {name_stream.str()};
using Mat_t = MaterialAnisotropic<dim, dim>;
using Sys_t = CellBase<dim, dim>;
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, std::vector<Real> input) -> Mat_t & {
return Mat_t::make(sys, n, input);
},
"cell"_a, "name"_a, "inputs"_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);
}
/**
* python binding for the optionally objective form of Anisotropic material
*/
template <Dim_t dim>
void add_material_orthotropic_helper(py::module & mod) {
std::stringstream name_stream{};
name_stream << "MaterialOrthotropic_" << dim << 'd';
const auto name {name_stream.str()};
using Mat_t = MaterialOrthotropic<dim, dim>;
using Sys_t = CellBase<dim, dim>;
using MatAniso_t = MaterialAnisotropic<dim,dim>;
// using MatBase_t = MaterialBase<dim,dim>;
py::class_<Mat_t, MatAniso_t>(mod, name.c_str())
.def_static("make",
[](Sys_t & sys, std::string n, std::vector<Real> input) -> Mat_t & {
return Mat_t::make(sys, n, input);
},
"cell"_a, "name"_a, "inputs"_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);
}
+
/**
* 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, MaterialBase<dim, dim>>(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, MaterialBase<dim, dim>>(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, MaterialBase<dim, dim>>(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, MaterialBase<dim, dim>>(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) {
std::stringstream name_stream{};
name_stream << "Material_" << dim << 'd';
const std::string name{name_stream.str()};
using Mat_t = MaterialBase<dim, dim>;
using FC_t = LocalFieldCollection<dim>;
using FCBase_t = FieldCollectionBase<dim, FC_t>;
py::class_<Mat_t>(mod, name.c_str()).
def_property_readonly
("collection",
[](Mat_t & material) -> FCBase_t &{
return material.get_collection();},
"returns the field collection containing internal "
"fields of this material",
py::return_value_policy::reference_internal);
- add_material_base_helper<dim>(mod);
+
+ //add_material_base_helper<dim>(mod);
add_material_anisotropic_helper<dim>(mod);
add_material_orthotropic_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_volume_calculator.hh b/src/common/intersection_volume_calculator.hh
index 7bc7737..080ec1d 100644
--- a/src/common/intersection_volume_calculator.hh
+++ b/src/common/intersection_volume_calculator.hh
@@ -1,207 +1,208 @@
/**
* @file intersection_volume_calculator.hh
*
* @author Ali Falsafi <ali.falsafi@epfl.ch>
*
* @date 04 June 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_VOLUME_CALCULATOR_H
#define INTERSECTION_VOLUME_CALCULATOR_H
//#gnclude <CGAL/Polyhedron_incremental_builder_3.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
#include <CGAL/IO/Polyhedron_iostream.h>
#include <CGAL/Nef_polyhedron_3.h>
#include <CGAL/IO/Nef_polyhedron_iostream_3.h>
#include <CGAL/Polyhedron_items_with_id_3.h>
#include <CGAL/Aff_transformation_3.h>
#include <CGAL/convex_hull_3.h>
#include <CGAL/Surface_mesh.h>
#include "common/ccoord_operations.hh"
#include "common/common.hh"
#include <vector>
#include <fstream>
#include <math.h>
namespace muSpectre {
typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
typedef CGAL::Polyhedron_3<Kernel, CGAL::Polyhedron_items_with_id_3> Polyhedron;
typedef CGAL::Nef_polyhedron_3<Kernel> Nef_polyhedron;
typedef Polyhedron::HalfedgeDS HalfedgeDS;
typedef Kernel::Point_3 Point_3;
typedef Kernel::Vector_3 Vector_3;
typedef Kernel::Tetrahedron_3 Tetrahedron;
typedef Polyhedron::Vertex_iterator Vertex_iterator;
typedef Polyhedron::Facet_iterator Facet_iterator;
typedef Polyhedron::Halfedge_around_facet_circulator Hafc;
typedef CGAL::Aff_transformation_3<Kernel> Transformation_mat;
using Dim_t = int;
using Real = double;
inline Polyhedron node_polyhedron_generator(std::array<Real, 3> origin, std::array<Real, 3> lengths) {
std::vector<Point_3> node_vertices ;
node_vertices.push_back(Point_3(origin[0] , origin[1] , origin[2] ));
node_vertices.push_back(Point_3(origin[0] + lengths[0] , origin[1] , origin[2] ));
node_vertices.push_back(Point_3(origin[0] , origin[1] + lengths[1], origin[2] ));
node_vertices.push_back(Point_3(origin[0] , origin[1] , origin[2] + lengths[2]));
node_vertices.push_back(Point_3(origin[0] , origin[1] + lengths[1], origin[2] + lengths[2]));
node_vertices.push_back(Point_3(origin[0] + lengths[0] , origin[1] , origin[2] + lengths[2]));
node_vertices.push_back(Point_3(origin[0] + lengths[0] , origin[1] + lengths[1], origin[2] ));
node_vertices.push_back(Point_3(origin[0] + lengths[0] , origin[1] + lengths[1], origin[2] + lengths[2]));
Polyhedron node_poly;
CGAL::convex_hull_3(node_vertices.begin(), node_vertices.end(), node_poly);
return node_poly;
}
inline Polyhedron convex_polyhedron_generator(std::vector<std::array<Real, 3>> convex_poly_vertices) {
std::vector<Point_3> poly_vertices ;
for (auto && point:convex_poly_vertices ){
poly_vertices.push_back(Point_3(point[0], point[1], point[2]));
}
Polyhedron convex_poly;
CGAL::convex_hull_3(poly_vertices.begin(), poly_vertices.end(), convex_poly);
return convex_poly;
}
template <Dim_t DimS>
class Correction{
public:
std::array<Real,3> correct_origin (std::array<Real,DimS> array);
std::array<Real,3> correct_length (std::array<Real,DimS> array );
std::vector<std::array<Real, 3>> correct_vector (std::vector<std::array<Real, DimS>> vector);
};
template<>
class Correction <3> {
public:
std::array<Real,3> correct_origin (std::array<Real,3> array){
return array;
}
std::array<Real,3> correct_length (std::array<Real,3> array){
return array;
}
std::vector<std::array<Real, 3>> correct_vector(std::vector<std::array<Real, 3>> vertices){
std::vector<std::array<Real, 3>> corrected_convex_poly_vertices;
return vertices;
}
};
template<>
class Correction <2> {
public:
std::vector<std::array<Real, 3>> correct_vector(std::vector<std::array<Real, 2>> vertices){
std::vector<std::array<Real, 3>> corrected_convex_poly_vertices;
for (auto && vertice : vertices){
- corrected_convex_poly_vertices.push_back({vertice.at(0), vertice.at(1), 0.0});
- corrected_convex_poly_vertices.push_back({vertice.at(0), vertice.at(1), 1.0});
+ corrected_convex_poly_vertices.push_back({vertice[0], vertice[1], 0.0});
+ corrected_convex_poly_vertices.push_back({vertice[0], vertice[1], 1.0});
}
return corrected_convex_poly_vertices;
}
std::array<Real,3> correct_origin (std::array<Real,2> array){
- return std::array<Real,3>{array.at(0),array.at(1),0.0};
+ return std::array<Real,3>{array[0],array[1],0.0};
}
std::array<Real,3> correct_length (std::array<Real,2> array){
- return std::array<Real,3>{array.at(0),array.at(1),1.0};
+ return std::array<Real,3>{array[0],array[1],1.0};
}
};
template <Dim_t DimS>
inline Real intersection_volume_calculator(std::vector<std::array<Real, DimS>> convex_poly_vertices,
- std::array<Real, DimS> origin, std::array<Real, DimS> lengths){
+ std::array<Real, DimS> origin,
+ std::array<Real, DimS> lengths){
Correction<DimS> correction;
std::vector<std::array<Real,3>> corrected_convex_poly_vertices (correction.correct_vector(convex_poly_vertices));
std::array<Real,3> corrected_origin(correction.correct_origin(origin));
std::array<Real,3> corrected_lengths(correction.correct_length(lengths));
std::vector<Point_3> result_vertex ;
std::vector<std::size_t> a_result_facet;
Vector_3 result_vector ;
- Point_3 center_point = Point_3( 0.0, 0.0, 0.0) ;
- Point_3 origin_coordinate_point (0,0,0), tet_vertices[4];
+ Point_3 center_point = Point_3(0.0, 0.0, 0.0) ;
+ Point_3 origin_coordinate_point (0, 0, 0), tet_vertices[4];
Real return_volume = 0.0;
Polyhedron node_poly, convex_poly, result_poly;
Tetrahedron tetra;
node_poly = node_polyhedron_generator(corrected_origin, corrected_lengths);
convex_poly = convex_polyhedron_generator(corrected_convex_poly_vertices);
Nef_polyhedron convex_poly_nef, node_poly_nef, result_nef ;
convex_poly_nef = Nef_polyhedron (convex_poly);
node_poly_nef = Nef_polyhedron (node_poly);
//computing the intersection:
result_nef = convex_poly_nef * node_poly_nef ;
if(result_nef.is_simple() and not result_nef.is_empty()) {
result_nef.convert_to_Polyhedron(result_poly);
//extracting vertex
int v_id = 0 ;
for ( Vertex_iterator v = result_poly.vertices_begin(); v != result_poly.vertices_end(); ++v){
result_vertex.push_back(v->point());
result_vector = v->point() - origin_coordinate_point;
center_point = center_point + result_vector;
v->id() = v_id++;
}
//calculating the center point of the resultant polyhedron:
Transformation_mat scale (1,0,0,0,0,1,0,0,0,0,1,0,v_id);
center_point = scale(center_point);
tet_vertices[0] = center_point;
int f_id = 0;
for ( Facet_iterator f = result_poly.facets_begin(); f != result_poly.facets_end(); ++f){
a_result_facet.clear();
f->id() = f_id++;
Hafc circ = f->facet_begin();
//assignig corresponding vertex to facets:
do {
a_result_facet.push_back(circ-> vertex()-> id());
} while ( ++circ != f->facet_begin());
//making tetreahedrals to calculate volumes:
for(int ii = 0 ; ii<3; ii++) tet_vertices[ii+1] = result_vertex[a_result_facet[ii]];
//calculating volumes:
tetra = Tetrahedron (tet_vertices[0], tet_vertices[1], tet_vertices[2], tet_vertices[3]);
return_volume += CGAL::to_double(tetra.volume());
}
}
return return_volume;
}
}
#endif //INTERSECTION_VOLUME_CALCULATOR
diff --git a/src/materials/material_linear_elastic1.hh b/src/materials/material_linear_elastic1.hh
index 35250ae..e15ead9 100644
--- a/src/materials/material_linear_elastic1.hh
+++ b/src/materials/material_linear_elastic1.hh
@@ -1,190 +1,189 @@
/**
* @file material_linear_elastic1.hh
*
* @author Till Junge <till.junge@epfl.ch>
*
* @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 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_LINEAR_ELASTIC1_H
#define MATERIAL_LINEAR_ELASTIC1_H
#include "common/common.hh"
#include "materials/material_muSpectre_base.hh"
#include "materials/materials_toolbox.hh"
namespace muSpectre {
template<Dim_t DimS, Dim_t DimM>
class MaterialLinearElastic1;
/**
* traits for objective linear elasticity
*/
template <Dim_t DimS, Dim_t DimM>
struct MaterialMuSpectre_traits<MaterialLinearElastic1<DimS, DimM>> {
using Parent = MaterialMuSpectre_traits<void>;//!< base for elasticity
//! global field collection
using GFieldCollection_t = typename
MaterialBase<DimS, DimM>::GFieldCollection_t;
//! expected map type for strain fields
using StrainMap_t = MatrixFieldMap<GFieldCollection_t, Real, DimM, DimM, true>;
//! expected map type for stress fields
using StressMap_t = MatrixFieldMap<GFieldCollection_t, Real, DimM, DimM>;
//! expected map type for tangent stiffness fields
using TangentMap_t = T4MatrixFieldMap<GFieldCollection_t, Real, DimM>;
//! 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<Dim_t DimS, Dim_t DimM>
class MaterialLinearElastic1:
public MaterialMuSpectre<MaterialLinearElastic1<DimS, DimM>, DimS, DimM>
{
public:
//! base class
using Parent = MaterialMuSpectre<MaterialLinearElastic1, DimS, DimM>;
/**
* 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 = T4Mat<Real, DimM>;
//! traits of this material
using traits = MaterialMuSpectre_traits<MaterialLinearElastic1>;
//! this law does not have any internal variables
using InternalVariables = typename traits::InternalVariables;
//! Hooke's law implementation
using Hooke = typename
MatTB::Hooke<DimM,
typename traits::StrainMap_t::reference,
typename traits::TangentMap_t::reference>;
//! 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 <class s_t>
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 <class s_t>
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 <Dim_t DimS, Dim_t DimM>
template <class s_t>
auto
MaterialLinearElastic1<DimS, DimM>::evaluate_stress(s_t && E)
-> decltype(auto) {
return Hooke::evaluate_stress(this->lambda, this->mu,
std::move(E));
}
/* ---------------------------------------------------------------------- */
template <Dim_t DimS, Dim_t DimM>
template <class s_t>
auto
MaterialLinearElastic1<DimS, DimM>::evaluate_stress_tangent(s_t && E)
-> decltype(auto) {
using Tangent_t = typename traits::TangentMap_t::reference;
-
// using mat = Eigen::Matrix<Real, DimM, DimM>;
// 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<double*>(this->C.data())),
// std::move(E)))} << std::endl;
return Hooke::evaluate_stress(this->lambda, this->mu,
Tangent_t(const_cast<double*>(this->C.data())),
std::move(E));
}
} // muSpectre
#endif /* MATERIAL_LINEAR_ELASTIC1_H */
diff --git a/src/materials/material_muSpectre_base.hh b/src/materials/material_muSpectre_base.hh
index 8a34ba8..fdbef3d 100644
--- a/src/materials/material_muSpectre_base.hh
+++ b/src/materials/material_muSpectre_base.hh
@@ -1,912 +1,912 @@
/**
* @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>
void add_pixel_split(const Ccoord_t<DimS> & pixel,
Real ratio,
InternalArgs ... args);
void add_pixel_split(const Ccoord_t<DimS> & pixel,
Real ratio = 1.0) override;
//add pixels intersecting to material to the material
void add_split_pixels_precipitate(std::vector<Ccoord_t<DimS>> intersected_pixles,
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 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_pixels,
std::vector<Real> intersection_ratios){
//assign precipitate materials:
for (auto && tup: akantu::zip(intersected_pixels, 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
+ // for the case that pixels are not split
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
+ // for the case that pixels are split
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 */