Diffusion µSpectre (feat/add_lam_pixels) Browse

add_lam_pixels

material_muSpectre_base.hh
View Options

			/**
			* @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 Lesser General Public License as
			* published by the Free Software Foundation, either version 3, or (at
			* your option) any later version.
			*
			* µSpectre is distributed in the hope that it will be useful, but
			* WITHOUT ANY WARRANTY; without even the implied warranty of
			* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
			* General Public License for more details.
			*
			* You should have received a copy of the GNU Lesser General Public License
			* along with µSpectre; see the file COPYING. If not, write to the
			* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
			* * Boston, MA 02111-1307, USA.
			*
			* Additional permission under GNU GPL version 3 section 7
			*
			* If you modify this Program, or any covered work, by linking or combining it
			* with proprietary FFT implementations or numerical libraries, containing parts
			* covered by the terms of those libraries' licenses, the licensors of this
			* Program grant you additional permission to convey the resulting work.
			*/

			#ifndef SRC_MATERIALS_MATERIAL_MUSPECTRE_BASE_HH_
			#define SRC_MATERIALS_MATERIAL_MUSPECTRE_BASE_HH_

			#include "common/common.hh"
			#include "materials/material_base.hh"
			#include "materials/materials_toolbox.hh"
			#include "laminate_homogenisation.hh"
			#include "common/field_collection.hh"
			#include "common/field.hh"
			#include "common//utilities.hh"
			#include "common/geometry.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>;
			// ! in order to save stress tensor of the pixels of a certian material
			using MStressField_t = typename Parent::MStressField_t;
			// ! in order to save stiffness tensor of the pixels of a certian material
			using MTangentField_t = typename Parent::MTangentField_t;

			//! Default constructor
			MaterialMuSpectre() = delete;

			//! Construct by name
			explicit MaterialMuSpectre(std::string name);

			//! Copy constructor
			MaterialMuSpectre(const MaterialMuSpectre & other) = delete;

			//! Move constructor
			MaterialMuSpectre(MaterialMuSpectre && other) = delete;

			//! Destructor
			virtual ~MaterialMuSpectre() = default;

			//! Factory
			template <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
			void initialise() override;

			// this fucntion is speicalized to assign partial material to a pixel
			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;
			//! computes stress
			void compute_stresses(const StrainField_t & F, StressField_t & P,
			Formulation form,
			SplitCell is_cell_split) override final;
			//! computes stress and tangent modulus
			void compute_stresses_tangent(const StrainField_t & F, StressField_t & P,
			TangentField_t & K, Formulation form,
			SplitCell is_cell_split) override final;

			protected:
			//! computes stress with the formulation available at compile time
			//! __attribute__ required by g++-6 and g++-7 because of this bug:
			//! https://gcc.gnu.org/bugzilla/show_bug.cgi?id=80947
			template <Formulation Form, SplitCell is_cell_split>
			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, SplitCell is_cell_split>
			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,
			SplitCell is_cell_split = SplitCell::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;
			mat_ref.allocate_optional_fields(cell.get_splitness());
			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.value().get().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.value().get().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,
			SplitCell is_cell_split) {
			switch (form) {
			case Formulation::finite_strain: {
			switch (is_cell_split) {
			case (SplitCell::no): {
			this->compute_stresses_worker<Formulation::finite_strain,
			SplitCell::no>(F, P);
			break;
			}
			case (SplitCell::simple): {
			this->compute_stresses_worker<Formulation::finite_strain,
			SplitCell::simple>(F, P);
			break;
			}
			case (SplitCell::laminate): {
			break;
			}
			}
			break;
			}
			case Formulation::small_strain: {
			switch (is_cell_split) {
			case (SplitCell::no): {
			this->compute_stresses_worker<Formulation::small_strain, SplitCell::no>(
			F, P);
			break;
			}
			case (SplitCell::simple): {
			this->compute_stresses_worker<Formulation::small_strain,
			SplitCell::simple>(F, P);
			break;
			}
			case (SplitCell::laminate): {
			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, SplitCell is_cell_split) {
			switch (form) {
			case Formulation::finite_strain: {
			switch (is_cell_split) {
			case (SplitCell::no): {
			this->compute_stresses_worker<Formulation::finite_strain,
			SplitCell::no>(F, P, K);
			break;
			}
			case (SplitCell::simple): {
			this->compute_stresses_worker<Formulation::finite_strain,
			SplitCell::simple>(F, P, K);
			break;
			}
			case (SplitCell::laminate): {
			break;
			}
			}
			break;
			}
			case Formulation::small_strain: {
			switch (is_cell_split) {
			case (SplitCell::no): {
			this->compute_stresses_worker<Formulation::small_strain, SplitCell::no>(
			F, P, K);
			break;
			}
			case (SplitCell::simple): {
			this->compute_stresses_worker<Formulation::small_strain,
			SplitCell::simple>(F, P, K);
			break;
			}
			case (SplitCell::laminate): {
			break;
			}
			}
			break;
			}
			default:
			throw std::runtime_error("Unknown formulation");
			break;
			}
			}

			/* ---------------------------------------------------------------------- */
			template <class Material, Dim_t DimS, Dim_t DimM>
			template <Formulation Form, SplitCell is_cell_split>
			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 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 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);
			};

			auto laminate_pixel = []() {};

			switch (is_cell_split) {
			case SplitCell::no: {
			non_split_pixel();
			break;
			}
			case SplitCell::simple: {
			split_pixel();
			break;
			}
			case SplitCell::laminate: {
			laminate_pixel();
			break;
			}
			}
			};

			auto constitutive_law_finite_strain = [this](Strains_t Strains,
			Stresses_t Stresses,
			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);
			};

			auto laminate_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)));
			*/
			};
			switch (is_cell_split) {
			case SplitCell::no: {
			non_split_pixel();
			break;
			}
			case SplitCell::simple: {
			split_pixel();
			break;
			}
			case SplitCell::laminate: {
			laminate_pixel();
			break;
			}
			}
			};

			iterable_proxy<NeedTangent::yes, is_cell_split> 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, SplitCell is_cell_split>
			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, 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 split 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);
			};

			auto lamiante_pixel = []() {};

			switch (is_cell_split) {
			case SplitCell::no: {
			non_split_pixel();
			break;
			}
			case SplitCell::simple: {
			split_pixel();
			break;
			}
			case SplitCell::laminate: {
			lamiante_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(junge): Figure this out: I can't std::move(internals...),
			// because if there are no internals, compilation fails with "no
			// matching function for call to ‘move()’'. These are tuples of
			// lvalue references, so it shouldn't be too bad, but still
			// irksome.

			// for the case that cell does not contain any split 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 split 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);
			};

			auto laminate_pixel = [&strain, &this_mat, ratio, &Stresses,
			internal_variables, &grad]() {};

			switch (is_cell_split) {
			case SplitCell::no: {
			non_split_pixel();
			break;
			}
			case SplitCell::simple: {
			split_pixel();
			break;
			}
			case SplitCell::laminate: {
			laminate_pixel();
			break;
			}
			}
			};

			iterable_proxy<NeedTangent::no, is_cell_split> 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, SplitCell is_cell_split>
			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, InternalReferences>;
			using iterator_category = std::forward_iterator_tag; //!< stl conformance

			//! Default constructor
			iterator() = delete;

			/** Iterator uses the material's internal variables field
			collection to iterate selectively over the global fields
			(such as the transformation gradient F and first
			Piola-Kirchhoff stress P.
			**/
			explicit iterator(const iterable_proxy & it, bool begin = true)
			: it{it}, strain_map{it.strain_field}, stress_map(it.stress_tup),
			index{begin ? 0 : it.material.internal_fields.size()} {}

			//! Copy constructor
			iterator(const iterator & other) = default;

			//! Move constructor
			iterator(iterator && other) = default;

			//! Destructor
			virtual ~iterator() = default;

			//! Copy assignment operator
			iterator & operator=(const iterator & other) = default;

			//! Move assignment operator
			iterator & operator=(iterator && other) = default;

			//! pre-increment
			inline iterator & operator++();
			//! dereference
			inline value_type operator*();
			//! inequality
			inline bool operator!=(const iterator & other) const;

			protected:
			const iterable_proxy & it; //!< ref to the proxy
			StrainMap_t strain_map; //!< map onto the global strain field
			//! map onto the global stress field and possibly tangent stiffness
			StressMapTup stress_map;
			size_t index; //!< index or pixel currently referred to

			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, SplitCell is_cell_split>
			bool
			MaterialMuSpectre<Material, DimS,
			DimM>::iterable_proxy<NeedTgt, is_cell_split>::iterator::
			operator!=(const iterator & other) const {
			return (this->index != other.index);
			}

			/* ---------------------------------------------------------------------- */
			template <class Material, Dim_t DimS, Dim_t DimM>
			template <MatTB::NeedTangent NeedTgt, SplitCell is_cell_split>
			typename MaterialMuSpectre<Material, DimS, DimM>::template iterable_proxy<
			NeedTgt, is_cell_split>::iterator &
			MaterialMuSpectre<Material, DimS,
			DimM>::iterable_proxy<NeedTgt, is_cell_split>::iterator::
			operator++() {
			this->index++;
			return *this;
			}

			/* ---------------------------------------------------------------------- */
			template <class Material, Dim_t DimS, Dim_t DimM>
			template <MatTB::NeedTangent NeedTgT, SplitCell is_cell_split>
			typename MaterialMuSpectre<Material, DimS, DimM>::template iterable_proxy<
			NeedTgT, is_cell_split>::iterator::value_type
			MaterialMuSpectre<Material, DimS, DimM>::iterable_proxy<
			NeedTgT, is_cell_split>::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 = 1.0;
			if (is_cell_split != SplitCell::no) {
			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));
			}
			} // namespace muSpectre

			#endif // SRC_MATERIALS_MATERIAL_MUSPECTRE_BASE_HH_

muSpectre/src/materials/material_muSpectre_base.hh096026929bfbfeat/add_lam_pixels

material_muSpectre_base.hhView Options

muSpectre/src/materials/material_muSpectre_base.hh
096026929bfbfeat/add_lam_pixels

material_muSpectre_base.hh
View Options