diff --git a/examples/io/dumper/locomotive_tools.hh b/examples/io/dumper/locomotive_tools.hh index e15f5d62b..cdf6d2c40 100644 --- a/examples/io/dumper/locomotive_tools.hh +++ b/examples/io/dumper/locomotive_tools.hh @@ -1,40 +1,40 @@ /** * @file locomotive_tools.hh * * @author Nicolas Richart * * @date creation: Fri Apr 13 2012 * @date last modification: Mon Aug 24 2015 * * @brief Header for the locomotive helper for the dumpers * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 of the License, or (at your option) any * later version. - * + * * Akantu 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 Lesser General Public License for more * details. - * + * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ void applyRotation(const ::akantu::Vector<::akantu::Real> & center, ::akantu::Real angle, const ::akantu::Array<::akantu::Real> & nodes, ::akantu::Array<::akantu::Real> & displacement, const ::akantu::Array<::akantu::UInt> & node_group); void fillColour(const ::akantu::Mesh & mesh, ::akantu::ElementTypeMapArray<::akantu::UInt> & colour); diff --git a/examples/new_material/local_material_damage.hh b/examples/new_material/local_material_damage.hh index 923e17a00..914d4f98e 100644 --- a/examples/new_material/local_material_damage.hh +++ b/examples/new_material/local_material_damage.hh @@ -1,119 +1,119 @@ /** * @file local_material_damage.hh * * @author Guillaume Anciaux * @author Marion Estelle Chambart * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Sep 29 2020 * * @brief Material isotropic elastic * * * @section LICENSE * * Copyright (©) 2015-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 of the License, or (at your option) any * later version. - * + * * Akantu 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 Lesser General Public License for more * details. - * + * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_LOCAL_MATERIAL_DAMAGE_HH_ #define AKANTU_LOCAL_MATERIAL_DAMAGE_HH_ namespace akantu { class LocalMaterialDamage : public Material { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: LocalMaterialDamage(SolidMechanicsModel & model, const ID & id = ""); virtual ~LocalMaterialDamage(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void initMaterial() override; /// constitutive law for all element of a type void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost) override; /// compute the potential energy for all elements void computePotentialEnergy(ElementType el_type) override; protected: /// constitutive law for a given quadrature point inline void computeStressOnQuad(Matrix & grad_u, Matrix & sigma, Real & damage); /// compute the potential energy for on element inline void computePotentialEnergyOnQuad(Matrix & grad_u, Matrix & sigma, Real & epot); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// compute the celerity of the fastest wave in the material inline Real getCelerity(const Element & element) const override; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Damage, damage, Real); private: /// the young modulus Real E; /// Poisson coefficient Real nu; /// First Lamé coefficient Real lambda; /// Second Lamé coefficient (shear modulus) Real mu; /// resistance to damage Real Yd; /// damage threshold Real Sd; /// Bulk modulus Real kpa; /// damage internal variable InternalField damage; }; } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "local_material_damage_inline_impl.hh" #endif /* AKANTU_LOCAL_MATERIAL_DAMAGE_HH_ */ diff --git a/examples/new_material/local_material_damage_inline_impl.hh b/examples/new_material/local_material_damage_inline_impl.hh index 646c7a9b9..aa5f42b0f 100644 --- a/examples/new_material/local_material_damage_inline_impl.hh +++ b/examples/new_material/local_material_damage_inline_impl.hh @@ -1,92 +1,92 @@ /** * @file local_material_damage_inline_impl.hh * * @author Guillaume Anciaux * @author Marion Estelle Chambart * @author Nicolas Richart * * @date creation: Wed Aug 04 2010 * @date last modification: Tue Sep 29 2020 * * @brief Implementation of the inline functions of the material damage * * * @section LICENSE * * Copyright (©) 2015-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 of the License, or (at your option) any * later version. - * + * * Akantu 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 Lesser General Public License for more * details. - * + * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef AKANTU_LOCAL_MATERIAL_DAMAGE_INLINE_IMPL_HH_ #define AKANTU_LOCAL_MATERIAL_DAMAGE_INLINE_IMPL_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ inline void LocalMaterialDamage::computeStressOnQuad(Matrix & grad_u, Matrix & sigma, Real & dam) { Real trace = grad_u.trace(); /// \sigma_{ij} = \lambda * (\nabla u)_{kk} * \delta_{ij} + \mu * (\nabla /// u_{ij} + \nabla u_{ji}) for (UInt i = 0; i < spatial_dimension; ++i) { for (UInt j = 0; j < spatial_dimension; ++j) { sigma(i, j) = (i == j) * lambda * trace + mu * (grad_u(i, j) + grad_u(j, i)); } } Real Y = 0; for (UInt i = 0; i < spatial_dimension; ++i) { for (UInt j = 0; j < spatial_dimension; ++j) { Y += sigma(i, j) * grad_u(i, j); } } Y *= 0.5; Real Fd = Y - Yd - Sd * dam; if (Fd > 0) dam = (Y - Yd) / Sd; dam = std::min(dam, 1.); sigma *= 1 - dam; } /* -------------------------------------------------------------------------- */ inline void LocalMaterialDamage::computePotentialEnergyOnQuad( Matrix & grad_u, Matrix & sigma, Real & epot) { epot = 0.; for (UInt i = 0, t = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j, ++t) epot += sigma(i, j) * (grad_u(i, j) - (i == j)); epot *= .5; } /* -------------------------------------------------------------------------- */ inline Real LocalMaterialDamage::getCelerity(__attribute__((unused)) const Element & element) const { // Here the fastest celerity is the push wave speed return (std::sqrt((2 * mu + lambda) / rho)); } } // namespace akantu #endif /* AKANTU_LOCAL_MATERIAL_DAMAGE_INLINE_IMPL_HH_ */ diff --git a/extra_packages/extra-materials/src/material_FE2/material_FE2.hh b/extra_packages/extra-materials/src/material_FE2/material_FE2.hh index 084f40da7..d886cfbc2 100644 --- a/extra_packages/extra-materials/src/material_FE2/material_FE2.hh +++ b/extra_packages/extra-materials/src/material_FE2/material_FE2.hh @@ -1,111 +1,110 @@ /** * @file material_FE2.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief Material for multi-scale simulations. It stores an * underlying RVE on each integration point of the material. * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "material_thermal.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MATERIAL_FE_2_HH_ #define AKANTU_MATERIAL_FE_2_HH_ namespace akantu { class SolidMechanicsModelRVE; } namespace akantu { /* -------------------------------------------------------------------------- */ /// /!\ This material works ONLY for meshes with a single element type!!!!! /* -------------------------------------------------------------------------- */ /** * MaterialFE2 * * parameters in the material files : * - mesh_file */ template class MaterialFE2 : public MaterialThermal { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ private: typedef MaterialThermal Parent; public: MaterialFE2(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialFE2(); typedef VoigtHelper voigt_h; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual void initMaterial(); /// constitutive law for all element of a type virtual void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); /// compute the tangent stiffness matrix for an element type - void computeTangentModuli(ElementType el_type, - Array & tangent_matrix, + void computeTangentModuli(ElementType el_type, Array & tangent_matrix, GhostType ghost_type = _not_ghost); /// advance alkali-silica reaction void advanceASR(const Matrix & prestrain); private: void initialize(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// Underlying RVE at each integration point std::vector> RVEs; /// Meshes for all RVEs std::vector> meshes; /// the element type of the associated mesh (this material handles only one /// type!!) ElementType el_type; /// the name of RVE mesh file ID mesh_file; /// Elastic stiffness tensor at each Gauss point (in voigt notation) InternalField C; /// number of gel pockets in each underlying RVE UInt nb_gel_pockets; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_FE2_inline_impl.hh" } // namespace akantu #endif /* AKANTU_MATERIAL_FE_2_HH_ */ diff --git a/extra_packages/extra-materials/src/material_damage/material_brittle.hh b/extra_packages/extra-materials/src/material_damage/material_brittle.hh index 5466e2c56..89df33608 100644 --- a/extra_packages/extra-materials/src/material_damage/material_brittle.hh +++ b/extra_packages/extra-materials/src/material_damage/material_brittle.hh @@ -1,110 +1,111 @@ /** * @file material_brittle.hh * * @author Aranda Ruiz Josue * @author Daniel Pino Muñoz * * * @brief Brittle damage law * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "material.hh" #include "material_damage.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MATERIAL_BRITTLE_HH_ #define AKANTU_MATERIAL_BRITTLE_HH_ namespace akantu { /** * Material brittle * * parameters in the material files : * - S_0 : Critical stress at low strain rate (default: 157e6) * - E_0 : Low strain rate threshold (default: 27e3) * - A,B,C,D : Fitting parameters for the critical stress at high strain * rates * (default: 1.622e-11, -1.3274e-6, 3.6544e-2, -181.38) */ template class MaterialBrittle : public MaterialDamage { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialBrittle(SolidMechanicsModel & model, const ID & id = ""); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void initMaterial() override; void updateInternalParameters() override; /// constitutive law for all element of a type - void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost) override; + void computeStress(ElementType el_type, + GhostType ghost_type = _not_ghost) override; protected: /// constitutive law for a given quadrature point inline void computeStressOnQuad(Matrix & grad_u, Matrix & grad_v, Matrix & sigma, Real & dam, Real & sigma_equivalent, Real & fracture_stress); inline void computeDamageAndStressOnQuad(Matrix & sigma, Real & dam, Real & sigma_c, Real & fracture_stress); /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ public: inline UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; inline void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; inline void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// strain rate arrays ordered by element types InternalField strain_rate_brittle; // polynome constants for critical stress value Real A; Real B; Real C; Real D; // minimum strain rate Real E_0; // Critical stress at low strain rates Real S_0; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_brittle_inline_impl.hh" } // namespace akantu #endif /* AKANTU_MATERIAL_brittle_HH_ */ diff --git a/extra_packages/extra-materials/src/material_damage/material_damage_iterative.hh b/extra_packages/extra-materials/src/material_damage/material_damage_iterative.hh index c844568f1..a9f2e05a7 100644 --- a/extra_packages/extra-materials/src/material_damage/material_damage_iterative.hh +++ b/extra_packages/extra-materials/src/material_damage/material_damage_iterative.hh @@ -1,145 +1,145 @@ /** * @file material_damage_iterative.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief Specialization of the class material damage to damage only one gauss * point at a time and propagate damage in a linear way. Max principal stress * criterion is used as a failure criterion. * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "material.hh" #include "material_damage.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MATERIAL_DAMAGE_ITERATIVE_HH_ #define AKANTU_MATERIAL_DAMAGE_ITERATIVE_HH_ namespace akantu { /** * Material damage iterative * * parameters in the material files : * - Sc */ template class MaterialDamageIterative : public MaterialDamage { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialDamageIterative(SolidMechanicsModel & model, const ID & id = ""); ~MaterialDamageIterative() override = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// virtual void updateInternalParameters(); void computeAllStresses(GhostType ghost_type = _not_ghost) override; /// update internal field damage virtual UInt updateDamage(); - UInt updateDamage(UInt quad_index, Real eq_stress, - ElementType el_type, GhostType ghost_type); + UInt updateDamage(UInt quad_index, Real eq_stress, ElementType el_type, + GhostType ghost_type); /// update energies after damage has been updated void updateEnergiesAfterDamage(ElementType el_type) override; /// compute the equivalent stress on each Gauss point (i.e. the max prinicpal /// stress) and normalize it by the tensile strength virtual void computeNormalizedEquivalentStress(const Array & grad_u, ElementType el_type, GhostType ghost_type = _not_ghost); /// find max normalized equivalent stress void findMaxNormalizedEquivalentStress(ElementType el_type, GhostType ghost_type = _not_ghost); protected: /// constitutive law for all element of a type void computeStress(ElementType el_type, - GhostType ghost_type = _not_ghost) override; + GhostType ghost_type = _not_ghost) override; inline void computeDamageAndStressOnQuad(Matrix & sigma, Real & dam); /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ inline UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; inline void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; inline void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get max normalized equivalent stress AKANTU_GET_MACRO(NormMaxEquivalentStress, norm_max_equivalent_stress, Real); /// get a non-const reference to the max normalized equivalent stress AKANTU_GET_MACRO_NOT_CONST(NormMaxEquivalentStress, norm_max_equivalent_stress, Real &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// resistance to damage RandomInternalField Sc; /// the reduction InternalField reduction_step; /// internal field to store equivalent stress on each Gauss point InternalField equivalent_stress; /// the number of total reductions steps until complete failure UInt max_reductions; /// damage increment Real prescribed_dam; /// maximum equivalent stress Real norm_max_equivalent_stress; /// deviation from max stress at which Gauss point will still get damaged Real dam_tolerance; /// define damage threshold at which damage will be set to 1 Real dam_threshold; /// maximum damage value Real max_damage; }; } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_damage_iterative_inline_impl.hh" #endif /* AKANTU_MATERIAL_DAMAGE_ITERATIVE_HH_ */ diff --git a/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage.hh b/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage.hh index 6690943d3..d15c70cbc 100644 --- a/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage.hh +++ b/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage.hh @@ -1,146 +1,145 @@ /** * @file material_orthotropic_damage.hh * @author Aurelia Cuba Ramos * @date Sun Mar 8 12:49:56 2015 * * @brief Material for orthotropic damage * * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 of the License, or (at your option) any * later version. * * Akantu 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "material_elastic.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MATERIAL_ORTHOTROPIC_DAMAGE_HH_ #define AKANTU_MATERIAL_ORTHOTROPIC_DAMAGE_HH_ namespace akantu { template class Parent = MaterialElastic> class MaterialOrthotropicDamage : public Parent { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialOrthotropicDamage(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialOrthotropicDamage(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void initMaterial() override; /// compute the tangent stiffness matrix for an element type - void computeTangentModuli(ElementType el_type, - Array & tangent_matrix, + void computeTangentModuli(ElementType el_type, Array & tangent_matrix, GhostType ghost_type = _not_ghost) override; protected: /// update the dissipated energy, must be called after the stress have been /// computed void updateEnergies(ElementType el_type) override; /// compute the tangent stiffness matrix for a given quadrature point inline void computeTangentModuliOnQuad( Matrix & tangent, const Matrix C, const Matrix & dam, const Matrix & dam_directions, Matrix & O_prime, Matrix & S_prime, Matrix & O, Matrix & S, Matrix & rotation_tmp); inline void computeDamageAndStressOnQuad(Matrix & sigma, Matrix & one_minus_D, Matrix & root_one_minus_D, Matrix & damage, Matrix & first_term, Matrix & third_term); /// rotate a Matrix of size dim*dim into the coordinate system of the FE /// computation inline void rotateIntoComputationFrame(const Matrix & to_rotate, Matrix & rotated, const Matrix & damage_directions, Matrix & rotation_tmp); /// rotate a Matrix of size dim*dim into the coordinate system of the damage inline void rotateIntoNewFrame(const Matrix & to_rotate, Matrix & rotated, const Matrix & damage_directions, Matrix & rotation_tmp); /// compute (1-D) inline void computeOneMinusD(Matrix & one_minus_D, const Matrix & damage); /// compute (1-D)^(1/2) inline void computeSqrtOneMinusD(const Matrix & one_minus_D, Matrix & sqrt_one_minus_D); /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// give the dissipated energy for the time step Real getDissipatedEnergy() const; // virtual Real getEnergy(std::string type); // virtual Real getEnergy(std::string energy_id, ElementType type, UInt index) // { // return Parent::getEnergy(energy_id, type, index); // }; AKANTU_GET_MACRO_NOT_CONST(Damage, damage, ElementTypeMapArray &); AKANTU_GET_MACRO(Damage, damage, const ElementTypeMapArray &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Damage, damage, Real) /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// damage internal variable InternalField damage; /// dissipated energy InternalField dissipated_energy; /// contain the current value of @f$ \int_0^{\epsilon}\sigma(\omega)d\omega /// @f$ the dissipated energy InternalField int_sigma; /// direction vectors for the damage frame InternalField damage_dir_vecs; Real eta; /// maximum damage value Real max_damage; }; } // namespace akantu #include "material_orthotropic_damage_tmpl.hh" #endif /* AKANTU_MATERIAL_ORTHOTROPIC_DAMAGE_HH_ */ diff --git a/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_tmpl.hh b/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_tmpl.hh index 2da4fcb7f..0d752a8d7 100644 --- a/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_tmpl.hh +++ b/extra_packages/extra-materials/src/material_damage/material_orthotropic_damage_tmpl.hh @@ -1,321 +1,320 @@ /** * @file material_orthotropic_damage_tmpl.hh * @author Aurelia Isabel Cuba Ramos * @date Sun Mar 8 12:54:30 2015 * * @brief Specialization of the material class for the orthotropic * damage material * * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 of the License, or (at your option) any * later version. * * Akantu 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_orthotropic_damage.hh" #include "solid_mechanics_model.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template class Parent> MaterialOrthotropicDamage::MaterialOrthotropicDamage( SolidMechanicsModel & model, const ID & id) : Parent(model, id), damage("damage", *this), dissipated_energy("damage dissipated energy", *this), int_sigma("integral of sigma", *this), damage_dir_vecs("damage_principal_directions", *this) { AKANTU_DEBUG_IN(); this->registerParam("eta", eta, 2., _pat_parsable | _pat_modifiable, "Damage sensitivity parameter"); this->registerParam("max_damage", max_damage, 0.99999, _pat_parsable | _pat_modifiable, "maximum damage value"); this->is_non_local = false; this->use_previous_stress = true; this->use_previous_gradu = true; /// use second order tensor for description of damage state this->damage.initialize(spatial_dimension * spatial_dimension); this->dissipated_energy.initialize(1); this->int_sigma.initialize(1); this->damage_dir_vecs.initialize(spatial_dimension * spatial_dimension); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class Parent> void MaterialOrthotropicDamage::initMaterial() { AKANTU_DEBUG_IN(); Parent::initMaterial(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the dissipated energy in each element by a trapezoidal approximation * of * @f$ Ed = \int_0^{\epsilon}\sigma(\omega)d\omega - * \frac{1}{2}\sigma:\epsilon@f$ */ template class Parent> void MaterialOrthotropicDamage::updateEnergies( ElementType el_type) { Parent::updateEnergies(el_type); this->computePotentialEnergy(el_type); auto epsilon_p = this->gradu.previous(el_type).begin(spatial_dimension, spatial_dimension); auto sigma_p = this->stress.previous(el_type).begin(spatial_dimension, spatial_dimension); auto epot = this->potential_energy(el_type).begin(); auto ints = this->int_sigma(el_type).begin(); auto ed = this->dissipated_energy(el_type).begin(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, _not_ghost); Matrix delta_gradu_it(grad_u); delta_gradu_it -= *epsilon_p; Matrix sigma_h(sigma); sigma_h += *sigma_p; Real dint = .5 * sigma_h.doubleDot(delta_gradu_it); *ints += dint; *ed = *ints - *epot; ++epsilon_p; ++sigma_p; ++epot; ++ints; ++ed; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; } /* -------------------------------------------------------------------------- */ template class Parent> void MaterialOrthotropicDamage::computeTangentModuli( - ElementType el_type, Array & tangent_matrix, - GhostType ghost_type) { + ElementType el_type, Array & tangent_matrix, GhostType ghost_type) { AKANTU_DEBUG_IN(); Parent::computeTangentModuli(el_type, tangent_matrix, ghost_type); /// get the damage array for current element type Array & dam = this->damage(el_type); auto dam_it = dam.begin(this->spatial_dimension, this->spatial_dimension); /// get the directions of damage for the current element type Array & dam_dirs = this->damage_dir_vecs(el_type); auto damage_directions_it = dam_dirs.begin(this->spatial_dimension, this->spatial_dimension); /// for the computation of the Cauchy stress the matrices (1-D) and /// (1-D)^(1/2) are needed. For the formulation see Engineering /// Damage Mechanics by Lemaitre and Desmorat. Matrix one_minus_D(this->spatial_dimension, this->spatial_dimension); Matrix sqrt_one_minus_D(this->spatial_dimension, this->spatial_dimension); Matrix one_minus_D_rot(spatial_dimension, spatial_dimension); Matrix sqrt_one_minus_D_rot(spatial_dimension, spatial_dimension); Matrix rotation_tmp(spatial_dimension, spatial_dimension); MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix); if (!(Math::are_float_equal((*dam_it).trace(), 0))) computeTangentModuliOnQuad(tangent, tangent, *dam_it, *damage_directions_it, one_minus_D, sqrt_one_minus_D, one_minus_D_rot, sqrt_one_minus_D_rot, rotation_tmp); ++dam_it; ++damage_directions_it; MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class Parent> void MaterialOrthotropicDamage:: computeTangentModuliOnQuad(Matrix & tangent, const Matrix C, const Matrix & dam, const Matrix & dam_directions, Matrix & O_prime, Matrix & S_prime, Matrix & O, Matrix & S, Matrix & rotation_tmp) { /// effect of damage on stiffness matrix: See Ragueneau et al. 2008, p. 423, /// ep. 7 Real trace_D = dam.trace(); this->computeOneMinusD(O_prime, dam); this->computeSqrtOneMinusD(O_prime, S_prime); this->rotateIntoComputationFrame(O_prime, O, dam_directions, rotation_tmp); this->rotateIntoComputationFrame(S_prime, S, dam_directions, rotation_tmp); /// compute new stiffness matrix in damage coordinate system if (spatial_dimension == 1) tangent *= (1 - dam(0, 0)); if (spatial_dimension == 2) { Real min_val = std::min((this->eta / spatial_dimension * trace_D), this->max_damage); /// first row tangent(0, 0) = (C(0, 0) * S(0, 0) * S(0, 0) + C(1, 0) * S(0, 1) * S(0, 1) - (min_val / 2. - 1. / 2) * (C(0, 0) + C(1, 0)) + (O(0, 0) * (C(0, 0) * O(0, 0) + C(1, 0) * O(1, 1))) / (trace_D - 2.)); tangent(0, 1) = (C(0, 1) * S(0, 0) * S(0, 0) + C(1, 1) * S(0, 1) * S(0, 1) - (min_val / 2. - 1. / 2) * (C(0, 1) + C(1, 1)) + (O(0, 0) * (C(0, 1) * O(0, 0) + C(1, 1) * O(1, 1))) / (trace_D - 2.)); tangent(0, 2) = (2. * C(2, 2) * S(0, 0) * S(0, 1) + (2. * C(2, 2) * O(0, 0) * O(0, 1)) / (trace_D - 2.)); /// second row tangent(1, 0) = (C(0, 0) * S(0, 1) * S(0, 1) + C(1, 0) * S(1, 1) * S(1, 1) - (min_val / 2. - 1. / 2) * (C(0, 0) + C(1, 0)) + (O(1, 1) * (C(0, 0) * O(0, 0) + C(1, 0) * O(1, 1))) / (trace_D - 2.)); tangent(1, 1) = (C(0, 1) * S(0, 1) * S(0, 1) + C(1, 1) * S(1, 1) * S(1, 1) - (min_val / 2. - 1. / 2) * (C(0, 1) + C(1, 1)) + (O(1, 1) * (C(0, 1) * O(0, 0) + C(1, 1) * O(1, 1))) / (trace_D - 2.)); tangent(1, 2) = (2. * C(2, 2) * S(0, 1) * S(1, 1) + (2. * C(2, 2) * O(0, 1) * O(1, 1)) / (trace_D - 2.)); /// third row tangent(2, 0) = ((O(0, 1) * (C(0, 0) * O(0, 0) + C(1, 0) * O(1, 1))) / (trace_D - 2.) + C(0, 0) * S(0, 0) * S(0, 1) + C(1, 0) * S(0, 1) * S(1, 1)); tangent(2, 1) = ((O(0, 1) * (C(0, 1) * O(0, 0) + C(1, 1) * O(1, 1))) / (trace_D - 2.) + C(0, 1) * S(0, 0) * S(0, 1) + C(1, 1) * S(0, 1) * S(1, 1)); tangent(2, 2) = ((2. * C(2, 2) * O(0, 1) * O(0, 1)) / (trace_D - 2.) + C(2, 2) * S(0, 1) * S(0, 1) + C(2, 2) * S(0, 0) * S(1, 1)); } if (spatial_dimension == 3) { } } /* -------------------------------------------------------------------------- */ template class Parent> inline void MaterialOrthotropicDamage:: computeDamageAndStressOnQuad(Matrix & sigma, Matrix & one_minus_D, Matrix & sqrt_one_minus_D, Matrix & damage, Matrix & first_term, Matrix & third_term) { /// Definition of Cauchy stress based on second order damage tensor: /// "Anisotropic damage modelling of biaxial behaviour and rupture /// of concrete strucutres", Ragueneau et al., 2008, Eq. 7 first_term.mul(sqrt_one_minus_D, sigma); first_term *= sqrt_one_minus_D; Real second_term = 0; for (UInt i = 0; i < this->spatial_dimension; ++i) { for (UInt j = 0; j < this->spatial_dimension; ++j) second_term += sigma(i, j) * one_minus_D(i, j); } second_term /= (this->spatial_dimension - damage.trace()); one_minus_D *= second_term; third_term.eye(1. / this->spatial_dimension * sigma.trace() * (1 - eta / (this->spatial_dimension) * damage.trace())); sigma.copy(first_term); sigma -= one_minus_D; sigma += third_term; } /* -------------------------------------------------------------------------- */ template class Parent> inline void MaterialOrthotropicDamage:: rotateIntoComputationFrame(const Matrix & to_rotate, Matrix & rotated, const Matrix & damage_directions, Matrix & rotation_tmp) { rotation_tmp.mul(to_rotate, damage_directions); rotated.mul(damage_directions, rotation_tmp); } /* -------------------------------------------------------------------------- */ template class Parent> inline void MaterialOrthotropicDamage::rotateIntoNewFrame( const Matrix & to_rotate, Matrix & rotated, const Matrix & damage_directions, Matrix & rotation_tmp) { rotation_tmp.mul(to_rotate, damage_directions); rotated.mul(damage_directions, rotation_tmp); } /* -------------------------------------------------------------------------- */ template class Parent> inline void MaterialOrthotropicDamage::computeOneMinusD( Matrix & one_minus_D, const Matrix & damage) { /// compute one minus one_minus_D.eye(); one_minus_D -= damage; } /* -------------------------------------------------------------------------- */ template class Parent> inline void MaterialOrthotropicDamage::computeSqrtOneMinusD( const Matrix & one_minus_D, Matrix & sqrt_one_minus_D) { /// To compute (1-D)^1/2 we need to check that we are in the /// principal coordinate system of the damage #ifndef AKANTU_NDEBUG for (UInt i = 0; i < this->spatial_dimension; ++i) { for (UInt j = 0; j < this->spatial_dimension; ++j) { if (i != j) AKANTU_DEBUG_ASSERT(Math::are_float_equal(one_minus_D(i, j), 0), "The damage tensor has off-diagonal parts"); } } #endif // AKANTU_NDEBUG /// compute (1-D)^1/2 sqrt_one_minus_D.copy(one_minus_D); for (UInt i = 0; i < this->spatial_dimension; ++i) sqrt_one_minus_D(i, i) = std::sqrt(sqrt_one_minus_D(i, i)); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/extra_packages/extra-materials/src/material_plastic/material_viscoplastic.hh b/extra_packages/extra-materials/src/material_plastic/material_viscoplastic.hh index eabfd8768..31976de37 100644 --- a/extra_packages/extra-materials/src/material_plastic/material_viscoplastic.hh +++ b/extra_packages/extra-materials/src/material_plastic/material_viscoplastic.hh @@ -1,99 +1,98 @@ /** * @file material_viscoplastic.hh * * @author Ramin Aghababaei * * * @brief Specialization of the material class for * MaterialLinearIsotropicHardening to include viscous effects (small * deformation) * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_voigthelper.hh" #include "material_plastic.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MATERIAL_VISCOPLASTIC_HH_ #define AKANTU_MATERIAL_VISCOPLASTIC_HH_ namespace akantu { /** * Material plastic isotropic * * parameters in the material files : * - h : Hardening parameter (default: 0) * - sigmay : Yield stress * - rate : Rate sensitivity * - edot0 : Reference strain rate * * - ts: Time step */ template class MaterialViscoPlastic : public MaterialPlastic { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialViscoPlastic(SolidMechanicsModel & model, const ID & id = ""); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// constitutive law for all element of a type virtual void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); /// compute the tangent stiffness matrix for an element type - void computeTangentModuli(ElementType el_type, - Array & tangent_matrix, + void computeTangentModuli(ElementType el_type, Array & tangent_matrix, GhostType ghost_type = _not_ghost); protected: inline void computeStressOnQuad(const Matrix & grad_u, const Matrix & previous_grad_u, Matrix & sigma, const Matrix & previous_sigma, Matrix & inelastic_strain, const Matrix & previous_inelastic_strain, Real & iso_hardening) const; inline void computeTangentModuliOnQuad( Matrix & tangent, const Matrix & grad_u, const Matrix & previous_grad_u, const Matrix & sigma_tensor, const Matrix & previous_sigma_tensor, const Real & iso_hardening) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// Rate sensitivity component (rate) Real rate; /// Reference strain rate (edot0) Real edot0; /// Time step (ts) Real ts; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_viscoplastic_inline_impl.hh" } // namespace akantu #endif /* AKANTU_MATERIAL_VISCOPLASTIC_HH_ */ diff --git a/extra_packages/extra-materials/src/material_viscoelastic/material_stiffness_proportional.hh b/extra_packages/extra-materials/src/material_viscoelastic/material_stiffness_proportional.hh index dbb3bac45..4b40c5f6f 100644 --- a/extra_packages/extra-materials/src/material_viscoelastic/material_stiffness_proportional.hh +++ b/extra_packages/extra-materials/src/material_viscoelastic/material_stiffness_proportional.hh @@ -1,98 +1,99 @@ /** * @file material_stiffness_proportional.hh * * @author David Simon Kammer * @author Nicolas Richart * * * @brief Material isotropic visco-elastic with viscosity proportional to the * stiffness * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "material.hh" #include "material_elastic.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MATERIAL_STIFFNESS_PROPORTIONAL_HH_ #define AKANTU_MATERIAL_STIFFNESS_PROPORTIONAL_HH_ namespace akantu { /** * Material visco-elastic @f[\sigma = E\epsilon + \alpha E* * \frac{d\epsilon}{dt}@f] * it can be seen as a Kelvin-Voigt solid with @f[\eta = \alpha E @f] * * The material satisfies the Caughey condition, the visco-elastic solid has the * same eigen-modes as the elastic one. (T.K. Caughey 1960 - Journal of Applied * Mechanics 27, 269-271. Classical normal modes in damped linear systems.) * * parameters in the material files : * - rho : density (default: 0) * - E : Young's modulus (default: 0) * - nu : Poisson's ratio (default: 1/2) * - Plane_Stress : if 0: plane strain, else: plane stress (default: 0) * - alpha : viscous ratio */ template class MaterialStiffnessProportional : public MaterialElastic { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialStiffnessProportional(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialStiffnessProportional(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void initMaterial() override; /// constitutive law for all element of a type - void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost) override; + void computeStress(ElementType el_type, + GhostType ghost_type = _not_ghost) override; /// compute the potential energy for all elements void computePotentialEnergy(ElementType el_type) override; protected: /// constitutive law for a given quadrature point // inline void computeStress(Real * F, Real * sigma); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// stress due to viscosity InternalField stress_viscosity; /// stress due to elasticity InternalField stress_elastic; /// viscous ratio Real alpha; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "material_elastic_caughey_inline_impl.hh" } // namespace akantu #endif /* AKANTU_MATERIAL_STIFFNESS_PROPORTIONAL_HH_ */ diff --git a/extra_packages/igfem/src/dumper_igfem_connectivity.hh b/extra_packages/igfem/src/dumper_igfem_connectivity.hh index 42191d61b..3f2c4c99e 100644 --- a/extra_packages/igfem/src/dumper_igfem_connectivity.hh +++ b/extra_packages/igfem/src/dumper_igfem_connectivity.hh @@ -1,116 +1,133 @@ /** * @file dumper_igfem_connectivity.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief Iterator for the IGFEM connectivity * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ #ifndef AKANTU_DUMPER_IGFEM_CONNECTIVITY_HH_ #define AKANTU_DUMPER_IGFEM_CONNECTIVITY_HH_ /* -------------------------------------------------------------------------- */ #include "dumper_igfem_element_iterator.hh" #include "dumper_igfem_generic_elemental_field.hh" /* -------------------------------------------------------------------------- */ namespace akantu { namespace dumpers { -/* -------------------------------------------------------------------------- */ - -template -class igfem_connectivity_field_iterator - : public igfem_element_iterator { - -public: - /* ------------------------------------------------------------------------ */ - /* Typedefs */ - /* ------------------------------------------------------------------------ */ - - typedef igfem_element_iterator - parent; - typedef typename types::return_type return_type; - typedef typename types::field_type field_type; - typedef typename types::array_iterator array_iterator; - -public: - /* ------------------------------------------------------------------------ */ - /* Constructors/Destructors */ - /* ------------------------------------------------------------------------ */ - - igfem_connectivity_field_iterator( - const field_type & field, const typename field_type::type_iterator & t_it, - const typename field_type::type_iterator & t_it_end, - const array_iterator & array_it, const array_iterator & array_it_end, - const GhostType ghost_type = _not_ghost, UInt sub_element = 0) - : parent(field, t_it, t_it_end, array_it, array_it_end, ghost_type, - sub_element) {} - - /* ------------------------------------------------------------------------ */ - /* Methods */ - /* ------------------------------------------------------------------------ */ - - return_type operator*() { - const Vector & element_connect = *this->array_it; - - /// get the local sub_element connectivity and the nodes per sub-element - UInt * sub_connec_ptr = - IGFEMHelper::getSubElementConnectivity(*this->tit, this->sub_element); - UInt nb_nodes_sub_el = - IGFEMHelper::getNbNodesPerSubElement(*this->tit, this->sub_element); - - /// get the global sub element connectivity - Vector sub_element_connect(nb_nodes_sub_el); - for (UInt i = 0; i < nb_nodes_sub_el; ++i) { - UInt lc = sub_connec_ptr[i]; - sub_element_connect(i) = element_connect(lc); + /* -------------------------------------------------------------------------- + */ + + template + class igfem_connectivity_field_iterator + : public igfem_element_iterator { + + public: + /* ------------------------------------------------------------------------ + */ + /* Typedefs */ + /* ------------------------------------------------------------------------ + */ + + typedef igfem_element_iterator + parent; + typedef typename types::return_type return_type; + typedef typename types::field_type field_type; + typedef typename types::array_iterator array_iterator; + + public: + /* ------------------------------------------------------------------------ + */ + /* Constructors/Destructors */ + /* ------------------------------------------------------------------------ + */ + + igfem_connectivity_field_iterator( + const field_type & field, + const typename field_type::type_iterator & t_it, + const typename field_type::type_iterator & t_it_end, + const array_iterator & array_it, const array_iterator & array_it_end, + const GhostType ghost_type = _not_ghost, UInt sub_element = 0) + : parent(field, t_it, t_it_end, array_it, array_it_end, ghost_type, + sub_element) {} + + /* ------------------------------------------------------------------------ + */ + /* Methods */ + /* ------------------------------------------------------------------------ + */ + + return_type operator*() { + const Vector & element_connect = *this->array_it; + + /// get the local sub_element connectivity and the nodes per sub-element + UInt * sub_connec_ptr = + IGFEMHelper::getSubElementConnectivity(*this->tit, this->sub_element); + UInt nb_nodes_sub_el = + IGFEMHelper::getNbNodesPerSubElement(*this->tit, this->sub_element); + + /// get the global sub element connectivity + Vector sub_element_connect(nb_nodes_sub_el); + for (UInt i = 0; i < nb_nodes_sub_el; ++i) { + UInt lc = sub_connec_ptr[i]; + sub_element_connect(i) = element_connect(lc); + } + + return sub_element_connect; } - return sub_element_connect; - } - - /* ------------------------------------------------------------------------ */ - /* Class Members */ - /* ------------------------------------------------------------------------ */ - -private: -}; - -/* -------------------------------------------------------------------------- */ -class IGFEMConnectivityField - : public IGFEMGenericElementalField, - igfem_connectivity_field_iterator> { - - /* ------------------------------------------------------------------------ */ - /* Typedefs */ - /* ------------------------------------------------------------------------ */ - -public: - typedef SingleType types; - typedef igfem_connectivity_field_iterator iterator; - typedef types::field_type field_type; - typedef IGFEMGenericElementalField - parent; - - /* ------------------------------------------------------------------------ */ - /* Constructors/Destructors */ - /* ------------------------------------------------------------------------ */ - - IGFEMConnectivityField(const field_type & field, - UInt spatial_dimension = _all_dimensions, - GhostType ghost_type = _not_ghost) - : parent(field, spatial_dimension, ghost_type) {} -}; - -/* -------------------------------------------------------------------------- */ + /* ------------------------------------------------------------------------ + */ + /* Class Members */ + /* ------------------------------------------------------------------------ + */ + + private: + }; + + /* -------------------------------------------------------------------------- + */ + class IGFEMConnectivityField + : public IGFEMGenericElementalField, + igfem_connectivity_field_iterator> { + + /* ------------------------------------------------------------------------ + */ + /* Typedefs */ + /* ------------------------------------------------------------------------ + */ + + public: + typedef SingleType types; + typedef igfem_connectivity_field_iterator iterator; + typedef types::field_type field_type; + typedef IGFEMGenericElementalField + parent; + + /* ------------------------------------------------------------------------ + */ + /* Constructors/Destructors */ + /* ------------------------------------------------------------------------ + */ + + IGFEMConnectivityField(const field_type & field, + UInt spatial_dimension = _all_dimensions, + GhostType ghost_type = _not_ghost) + : parent(field, spatial_dimension, ghost_type) {} + }; + + /* -------------------------------------------------------------------------- + */ } // namespace dumpers } // namespace akantu /* -------------------------------------------------------------------------- */ #endif /*AKANTU_DUMPER_IGFEM_CONNECTIVITY_HH_ */ diff --git a/extra_packages/igfem/src/dumper_igfem_element_iterator.hh b/extra_packages/igfem/src/dumper_igfem_element_iterator.hh index 09a320438..a044f8522 100644 --- a/extra_packages/igfem/src/dumper_igfem_element_iterator.hh +++ b/extra_packages/igfem/src/dumper_igfem_element_iterator.hh @@ -1,181 +1,199 @@ /** * @file dumper_igfem_element_iterator.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief Iterators for IGFEM elemental fields * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ #ifndef AKANTU_DUMPER_IGFEM_ELEMENT_ITERATOR_HH_ #define AKANTU_DUMPER_IGFEM_ELEMENT_ITERATOR_HH_ /* -------------------------------------------------------------------------- */ #include "element.hh" #include "igfem_helper.hh" /* -------------------------------------------------------------------------- */ namespace akantu { namespace dumpers { -/* -------------------------------------------------------------------------- */ + /* -------------------------------------------------------------------------- + */ + + template class final_iterator> + class igfem_element_iterator { + /* ------------------------------------------------------------------------ + */ + /* Typedefs */ + /* ------------------------------------------------------------------------ + */ + public: + typedef typename types::it_type it_type; + typedef typename types::field_type field_type; + typedef typename types::array_type array_type; + typedef typename types::array_iterator array_iterator; + typedef final_iterator iterator; + + public: + /* ------------------------------------------------------------------------ + */ + /* Constructors/Destructors */ + /* ------------------------------------------------------------------------ + */ + + igfem_element_iterator(const field_type & field, + const typename field_type::type_iterator & t_it, + const typename field_type::type_iterator & t_it_end, + const array_iterator & array_it, + const array_iterator & array_it_end, + const GhostType ghost_type = _not_ghost, + UInt sub_element = 0) + : field(field), tit(t_it), tit_end(t_it_end), array_it(array_it), + array_it_end(array_it_end), ghost_type(ghost_type), + sub_element(sub_element) {} + + /* ------------------------------------------------------------------------ + */ + /* Methods */ + /* ------------------------------------------------------------------------ + */ + + public: + bool operator!=(const iterator & it) const { + return (ghost_type != it.ghost_type) || + (tit != it.tit || + ((array_it != it.array_it) || sub_element != it.sub_element)); + } -template class final_iterator> -class igfem_element_iterator { - /* ------------------------------------------------------------------------ */ - /* Typedefs */ - /* ------------------------------------------------------------------------ */ -public: - typedef typename types::it_type it_type; - typedef typename types::field_type field_type; - typedef typename types::array_type array_type; - typedef typename types::array_iterator array_iterator; - typedef final_iterator iterator; - -public: - /* ------------------------------------------------------------------------ */ - /* Constructors/Destructors */ - /* ------------------------------------------------------------------------ */ - - igfem_element_iterator(const field_type & field, - const typename field_type::type_iterator & t_it, - const typename field_type::type_iterator & t_it_end, - const array_iterator & array_it, - const array_iterator & array_it_end, - const GhostType ghost_type = _not_ghost, - UInt sub_element = 0) - : field(field), tit(t_it), tit_end(t_it_end), array_it(array_it), - array_it_end(array_it_end), ghost_type(ghost_type), - sub_element(sub_element) {} - - /* ------------------------------------------------------------------------ */ - /* Methods */ - /* ------------------------------------------------------------------------ */ - -public: - bool operator!=(const iterator & it) const { - return (ghost_type != it.ghost_type) || - (tit != it.tit || - ((array_it != it.array_it) || sub_element != it.sub_element)); - } - - iterator & operator++() { - if (!this->sub_element) - this->sub_element += 1; - else { - ++array_it; - this->sub_element = 0; - while (array_it == array_it_end && tit != tit_end) { - ++tit; - if (tit != tit_end) { - const array_type & vect = field(*tit, ghost_type); - UInt _nb_data_per_elem = getNbDataPerElem(*tit); - UInt nb_component = vect.getNbComponent(); - UInt size = (vect.getSize() * nb_component) / _nb_data_per_elem; - - array_it = vect.begin_reinterpret(_nb_data_per_elem, size); - array_it_end = vect.end_reinterpret(_nb_data_per_elem, size); + iterator & operator++() { + if (!this->sub_element) + this->sub_element += 1; + else { + ++array_it; + this->sub_element = 0; + while (array_it == array_it_end && tit != tit_end) { + ++tit; + if (tit != tit_end) { + const array_type & vect = field(*tit, ghost_type); + UInt _nb_data_per_elem = getNbDataPerElem(*tit); + UInt nb_component = vect.getNbComponent(); + UInt size = (vect.getSize() * nb_component) / _nb_data_per_elem; + + array_it = vect.begin_reinterpret(_nb_data_per_elem, size); + array_it_end = vect.end_reinterpret(_nb_data_per_elem, size); + } } } + return *(static_cast(this)); } - return *(static_cast(this)); - } - - ElementType getType() { - ElementType sub_type = IGFEMHelper::getSubElementType(*tit, sub_element); - return sub_type; - } - - /// get IOHelperType for sub-element - UInt element_type() { return getIOHelperType(this->getType()); } - - /// get current parent element???? - Element getCurrentElement() { - return Element(*tit, array_it.getCurrentIndex()); - } - - UInt getNbDataPerElem(ElementType type) const { - /// nb of data per parent element! - if (!nb_data_per_elem.exists(type, ghost_type)) - return field(type, ghost_type).getNbComponent(); - - return nb_data_per_elem(type, ghost_type); - } - - void setNbDataPerElem(const ElementTypeMap & nb_data) { - /// nb of data per parent element! - this->nb_data_per_elem = nb_data; - } - - /* ------------------------------------------------------------------------ */ - /* Class Members */ - /* ------------------------------------------------------------------------ */ - -protected: - /// the field to iterate on - const field_type & field; - /// field iterator - typename field_type::type_iterator tit; - /// field iterator end - typename field_type::type_iterator tit_end; - /// array iterator - array_iterator array_it; - /// internal iterator end - array_iterator array_it_end; - /// ghost type identification - const GhostType ghost_type; - /// number of data per element - ElementTypeMap nb_data_per_elem; - /// index of sub-element - UInt sub_element; - /// sub_element end - UInt sub_element_end; -}; -/* -------------------------------------------------------------------------- */ -template -class igfem_elemental_field_iterator - : public igfem_element_iterator { -public: - /* ------------------------------------------------------------------------ */ - /* Typedefs */ - /* ------------------------------------------------------------------------ */ - - typedef igfem_element_iterator< - types, ::akantu::dumpers::igfem_elemental_field_iterator> - parent; - typedef typename types::it_type it_type; - typedef typename types::return_type return_type; - typedef typename types::field_type field_type; - typedef typename types::array_iterator array_iterator; - -public: - /* ------------------------------------------------------------------------ */ - /* Constructors/Destructors */ - /* ------------------------------------------------------------------------ */ - - igfem_elemental_field_iterator( - const field_type & field, const typename field_type::type_iterator & t_it, - const typename field_type::type_iterator & t_it_end, - const array_iterator & array_it, const array_iterator & array_it_end, - const GhostType ghost_type = _not_ghost, UInt sub_element = 0) - : parent(field, t_it, t_it_end, array_it, array_it_end, ghost_type, - sub_element) {} - - /* ------------------------------------------------------------------------ */ - /* Methods */ - /* ------------------------------------------------------------------------ */ - - return_type operator*() { return *this->array_it; } - -private: -}; + ElementType getType() { + ElementType sub_type = IGFEMHelper::getSubElementType(*tit, sub_element); + return sub_type; + } -/* -------------------------------------------------------------------------- */ + /// get IOHelperType for sub-element + UInt element_type() { return getIOHelperType(this->getType()); } + + /// get current parent element???? + Element getCurrentElement() { + return Element(*tit, array_it.getCurrentIndex()); + } + + UInt getNbDataPerElem(ElementType type) const { + /// nb of data per parent element! + if (!nb_data_per_elem.exists(type, ghost_type)) + return field(type, ghost_type).getNbComponent(); + + return nb_data_per_elem(type, ghost_type); + } + + void setNbDataPerElem(const ElementTypeMap & nb_data) { + /// nb of data per parent element! + this->nb_data_per_elem = nb_data; + } + + /* ------------------------------------------------------------------------ + */ + /* Class Members */ + /* ------------------------------------------------------------------------ + */ + + protected: + /// the field to iterate on + const field_type & field; + /// field iterator + typename field_type::type_iterator tit; + /// field iterator end + typename field_type::type_iterator tit_end; + /// array iterator + array_iterator array_it; + /// internal iterator end + array_iterator array_it_end; + /// ghost type identification + const GhostType ghost_type; + /// number of data per element + ElementTypeMap nb_data_per_elem; + /// index of sub-element + UInt sub_element; + /// sub_element end + UInt sub_element_end; + }; + + /* -------------------------------------------------------------------------- + */ + template + class igfem_elemental_field_iterator + : public igfem_element_iterator { + public: + /* ------------------------------------------------------------------------ + */ + /* Typedefs */ + /* ------------------------------------------------------------------------ + */ + + typedef igfem_element_iterator< + types, ::akantu::dumpers::igfem_elemental_field_iterator> + parent; + typedef typename types::it_type it_type; + typedef typename types::return_type return_type; + typedef typename types::field_type field_type; + typedef typename types::array_iterator array_iterator; + + public: + /* ------------------------------------------------------------------------ + */ + /* Constructors/Destructors */ + /* ------------------------------------------------------------------------ + */ + + igfem_elemental_field_iterator( + const field_type & field, + const typename field_type::type_iterator & t_it, + const typename field_type::type_iterator & t_it_end, + const array_iterator & array_it, const array_iterator & array_it_end, + const GhostType ghost_type = _not_ghost, UInt sub_element = 0) + : parent(field, t_it, t_it_end, array_it, array_it_end, ghost_type, + sub_element) {} + + /* ------------------------------------------------------------------------ + */ + /* Methods */ + /* ------------------------------------------------------------------------ + */ + + return_type operator*() { return *this->array_it; } + + private: + }; + + /* -------------------------------------------------------------------------- + */ } // namespace dumpers } // namespace akantu /* -------------------------------------------------------------------------- */ #endif /* AKANTU_DUMPER_IGFEM_ELEMENT_ITERATOR_HH_ */ diff --git a/extra_packages/igfem/src/dumper_igfem_element_partition.hh b/extra_packages/igfem/src/dumper_igfem_element_partition.hh index 88acea1ef..e098b69d0 100644 --- a/extra_packages/igfem/src/dumper_igfem_element_partition.hh +++ b/extra_packages/igfem/src/dumper_igfem_element_partition.hh @@ -1,103 +1,122 @@ /** * @file dumper_igfem_element_partition.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief Element partition field for IGFEM sub-elements * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ namespace akantu { namespace dumpers { -/* -------------------------------------------------------------------------- */ -template -class igfem_element_partition_field_iterator - : public igfem_element_iterator { + /* -------------------------------------------------------------------------- + */ + template + class igfem_element_partition_field_iterator + : public igfem_element_iterator { - /* ------------------------------------------------------------------------ */ - /* Typedefs */ - /* ------------------------------------------------------------------------ */ -public: - typedef igfem_element_iterator - parent; - typedef typename types::return_type return_type; - typedef typename types::array_iterator array_iterator; - typedef typename types::field_type field_type; + /* ------------------------------------------------------------------------ + */ + /* Typedefs */ + /* ------------------------------------------------------------------------ + */ + public: + typedef igfem_element_iterator< + types, dumpers::igfem_element_partition_field_iterator> + parent; + typedef typename types::return_type return_type; + typedef typename types::array_iterator array_iterator; + typedef typename types::field_type field_type; - /* ------------------------------------------------------------------------ */ - /* Constructors/Destructors */ - /* ------------------------------------------------------------------------ */ -public: - igfem_element_partition_field_iterator( - const field_type & field, const typename field_type::type_iterator & t_it, - const typename field_type::type_iterator & t_it_end, - const array_iterator & array_it, const array_iterator & array_it_end, - const GhostType ghost_type = _not_ghost, UInt sub_element = 0) - : parent(field, t_it, t_it_end, array_it, array_it_end, ghost_type, - sub_element) { - prank = StaticCommunicator::getStaticCommunicator().whoAmI(); - } + /* ------------------------------------------------------------------------ + */ + /* Constructors/Destructors */ + /* ------------------------------------------------------------------------ + */ + public: + igfem_element_partition_field_iterator( + const field_type & field, + const typename field_type::type_iterator & t_it, + const typename field_type::type_iterator & t_it_end, + const array_iterator & array_it, const array_iterator & array_it_end, + const GhostType ghost_type = _not_ghost, UInt sub_element = 0) + : parent(field, t_it, t_it_end, array_it, array_it_end, ghost_type, + sub_element) { + prank = StaticCommunicator::getStaticCommunicator().whoAmI(); + } - /* ------------------------------------------------------------------------ */ - /* Methods */ - /* ------------------------------------------------------------------------ */ -public: - return_type operator*() { return return_type(1, prank); } + /* ------------------------------------------------------------------------ + */ + /* Methods */ + /* ------------------------------------------------------------------------ + */ + public: + return_type operator*() { return return_type(1, prank); } - /* ------------------------------------------------------------------------ */ - /* Class Members */ - /* ------------------------------------------------------------------------ */ -protected: - UInt prank; -}; + /* ------------------------------------------------------------------------ + */ + /* Class Members */ + /* ------------------------------------------------------------------------ + */ + protected: + UInt prank; + }; -/* -------------------------------------------------------------------------- */ -template -class IGFEMElementPartitionField : public IGFEMGenericElementalField< - SingleType, - igfem_element_partition_field_iterator> { -public: - /* ------------------------------------------------------------------------ */ - /* Typedefs */ - /* ------------------------------------------------------------------------ */ + /* -------------------------------------------------------------------------- + */ + template + class IGFEMElementPartitionField + : public IGFEMGenericElementalField< + SingleType, + igfem_element_partition_field_iterator> { + public: + /* ------------------------------------------------------------------------ + */ + /* Typedefs */ + /* ------------------------------------------------------------------------ + */ - typedef SingleType types; - typedef igfem_element_partition_field_iterator iterator; - typedef IGFEMGenericElementalField - parent; - typedef typename types::field_type field_type; + typedef SingleType types; + typedef igfem_element_partition_field_iterator iterator; + typedef IGFEMGenericElementalField + parent; + typedef typename types::field_type field_type; -public: - /* ------------------------------------------------------------------------ */ - /* Constructors/Destructors */ - /* ------------------------------------------------------------------------ */ + public: + /* ------------------------------------------------------------------------ + */ + /* Constructors/Destructors */ + /* ------------------------------------------------------------------------ + */ - IGFEMElementPartitionField(const field_type & field, - UInt spatial_dimension = _all_dimensions, - GhostType ghost_type = _not_ghost, - ElementKind kind = _ek_igfem) - : parent(field, spatial_dimension, ghost_type, kind) { - this->homogeneous = true; - } + IGFEMElementPartitionField(const field_type & field, + UInt spatial_dimension = _all_dimensions, + GhostType ghost_type = _not_ghost, + ElementKind kind = _ek_igfem) + : parent(field, spatial_dimension, ghost_type, kind) { + this->homogeneous = true; + } - /* ------------------------------------------------------------------------ */ - /* Methods */ - /* ------------------------------------------------------------------------ */ + /* ------------------------------------------------------------------------ + */ + /* Methods */ + /* ------------------------------------------------------------------------ + */ - UInt getDim() { return 1; } -}; + UInt getDim() { return 1; } + }; -/* -------------------------------------------------------------------------- */ + /* -------------------------------------------------------------------------- + */ } // namespace dumpers } // namespace akantu diff --git a/extra_packages/igfem/src/dumper_igfem_elemental_field.hh b/extra_packages/igfem/src/dumper_igfem_elemental_field.hh index 323c2dc7e..c3e53d4f6 100644 --- a/extra_packages/igfem/src/dumper_igfem_elemental_field.hh +++ b/extra_packages/igfem/src/dumper_igfem_elemental_field.hh @@ -1,56 +1,61 @@ /** * @file dumper_igfem_elemental_field.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief description of IGFEM elemental fields * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ #ifndef AKANTU_DUMPER_IGFEM_ELEMENTAL_FIELD_HH_ #define AKANTU_DUMPER_IGFEM_ELEMENTAL_FIELD_HH_ /* -------------------------------------------------------------------------- */ #include "dumper_field.hh" #include "dumper_igfem_generic_elemental_field.hh" #include "static_communicator.hh" /* -------------------------------------------------------------------------- */ namespace akantu { namespace dumpers { -/* -------------------------------------------------------------------------- */ - -template class ret = Vector, - bool filtered = false> -class IGFEMElementalField - : public IGFEMGenericElementalField, - igfem_elemental_field_iterator> { - -public: - /* ------------------------------------------------------------------------ */ - /* Typedefs */ - /* ------------------------------------------------------------------------ */ - - typedef SingleType types; - typedef typename types::field_type field_type; - typedef elemental_field_iterator iterator; - - /* ------------------------------------------------------------------------ */ - /* Constructors/Destructors */ - /* ------------------------------------------------------------------------ */ - - IGFEMElementalField(const field_type & field, - UInt spatial_dimension = _all_dimensions, - GhostType ghost_type = _not_ghost, - ElementKind element_kind = _ek_igfem) - : IGFEMGenericElementalField( - field, spatial_dimension, ghost_type, element_kind) {} -}; + /* -------------------------------------------------------------------------- + */ + + template class ret = Vector, + bool filtered = false> + class IGFEMElementalField + : public IGFEMGenericElementalField, + igfem_elemental_field_iterator> { + + public: + /* ------------------------------------------------------------------------ + */ + /* Typedefs */ + /* ------------------------------------------------------------------------ + */ + + typedef SingleType types; + typedef typename types::field_type field_type; + typedef elemental_field_iterator iterator; + + /* ------------------------------------------------------------------------ + */ + /* Constructors/Destructors */ + /* ------------------------------------------------------------------------ + */ + + IGFEMElementalField(const field_type & field, + UInt spatial_dimension = _all_dimensions, + GhostType ghost_type = _not_ghost, + ElementKind element_kind = _ek_igfem) + : IGFEMGenericElementalField( + field, spatial_dimension, ghost_type, element_kind) {} + }; } // namespace dumpers } // namespace akantu #endif /* AKANTU_DUMPER_IGFEM_ELEMENTAL_FIELD_HH_ */ diff --git a/extra_packages/igfem/src/dumper_igfem_generic_elemental_field.hh b/extra_packages/igfem/src/dumper_igfem_generic_elemental_field.hh index 66bb3bf07..b88a98a72 100644 --- a/extra_packages/igfem/src/dumper_igfem_generic_elemental_field.hh +++ b/extra_packages/igfem/src/dumper_igfem_generic_elemental_field.hh @@ -1,143 +1,152 @@ /** * @file dumper_igfem_generic_elemental_field.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief generic interface IGFEM elemental fields * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #ifndef AKANTU_DUMPER_IGFEM_GENERIC_ELEMENTAL_FIELD_HH_ #define AKANTU_DUMPER_IGFEM_GENERIC_ELEMENTAL_FIELD_HH_ /* -------------------------------------------------------------------------- */ #include "dumper_generic_elemental_field.hh" #include "dumper_igfem_element_iterator.hh" /* -------------------------------------------------------------------------- */ namespace akantu { namespace dumpers { -/* -------------------------------------------------------------------------- */ -template class iterator_type> -class IGFEMGenericElementalField - : public GenericElementalField<_types, iterator_type> { - - /* ------------------------------------------------------------------------ */ - /* Typedefs */ - /* ------------------------------------------------------------------------ */ - -public: - typedef _types types; - typedef typename types::data_type data_type; - typedef typename types::it_type it_type; - typedef typename types::field_type field_type; - typedef typename types::array_type array_type; - typedef typename types::array_iterator array_iterator; - typedef typename field_type::type_iterator field_type_iterator; - typedef iterator_type iterator; - - /* ------------------------------------------------------------------------ */ - /* Constructors/Destructors */ - /* ------------------------------------------------------------------------ */ - -public: - IGFEMGenericElementalField(const field_type & field, - UInt spatial_dimension = _all_dimensions, - GhostType ghost_type = _not_ghost, - ElementKind kind = _ek_igfem) - : - - GenericElementalField(field, spatial_dimension, - ghost_type, kind) { - this->checkHomogeneity(); - } - - /* ------------------------------------------------------------------------ */ - /* Methods */ - /* ------------------------------------------------------------------------ */ - -public: - /// return the size of the contained data: i.e. the number of elements ? - virtual UInt size() { - this->checkHomogeneity(); - return ((this->nb_total_element) * 2); - } - - virtual iterator begin() { - field_type_iterator tit; - field_type_iterator end; - UInt sub_element = 0; - - /// type iterators on the elemental field - tit = this->field.firstType(this->spatial_dimension, this->ghost_type, - this->element_kind); - end = this->field.lastType(this->spatial_dimension, this->ghost_type, - this->element_kind); - - /// skip all types without data - ElementType type = *tit; - for (; tit != end && this->field(*tit, this->ghost_type).getSize() == 0; - ++tit) { + /* -------------------------------------------------------------------------- + */ + template class iterator_type> + class IGFEMGenericElementalField + : public GenericElementalField<_types, iterator_type> { + + /* ------------------------------------------------------------------------ + */ + /* Typedefs */ + /* ------------------------------------------------------------------------ + */ + + public: + typedef _types types; + typedef typename types::data_type data_type; + typedef typename types::it_type it_type; + typedef typename types::field_type field_type; + typedef typename types::array_type array_type; + typedef typename types::array_iterator array_iterator; + typedef typename field_type::type_iterator field_type_iterator; + typedef iterator_type iterator; + + /* ------------------------------------------------------------------------ + */ + /* Constructors/Destructors */ + /* ------------------------------------------------------------------------ + */ + + public: + IGFEMGenericElementalField(const field_type & field, + UInt spatial_dimension = _all_dimensions, + GhostType ghost_type = _not_ghost, + ElementKind kind = _ek_igfem) + : + + GenericElementalField(field, spatial_dimension, + ghost_type, kind) { + this->checkHomogeneity(); + } + + /* ------------------------------------------------------------------------ + */ + /* Methods */ + /* ------------------------------------------------------------------------ + */ + + public: + /// return the size of the contained data: i.e. the number of elements ? + virtual UInt size() { + this->checkHomogeneity(); + return ((this->nb_total_element) * 2); } - type = *tit; - - if (tit == end) - return this->end(); - - /// getting information for the field of the given type - const array_type & vect = this->field(type, this->ghost_type); - UInt nb_data_per_elem = this->getNbDataPerElem(type); - UInt nb_component = vect.getNbComponent(); - UInt size = (vect.getSize() * nb_component) / nb_data_per_elem; - - /// define element-wise iterator - array_iterator it = vect.begin_reinterpret(nb_data_per_elem, size); - array_iterator it_end = vect.end_reinterpret(nb_data_per_elem, size); - /// define data iterator - iterator rit = iterator(this->field, tit, end, it, it_end, this->ghost_type, - sub_element); - rit.setNbDataPerElem(this->nb_data_per_elem); - return rit; - } - - virtual iterator end() { - field_type_iterator tit; - field_type_iterator end; - UInt sub_element = 0; - - tit = this->field.firstType(this->spatial_dimension, this->ghost_type, - this->element_kind); - end = this->field.lastType(this->spatial_dimension, this->ghost_type, - this->element_kind); - - ElementType type = *tit; - for (; tit != end; ++tit) + + virtual iterator begin() { + field_type_iterator tit; + field_type_iterator end; + UInt sub_element = 0; + + /// type iterators on the elemental field + tit = this->field.firstType(this->spatial_dimension, this->ghost_type, + this->element_kind); + end = this->field.lastType(this->spatial_dimension, this->ghost_type, + this->element_kind); + + /// skip all types without data + ElementType type = *tit; + for (; tit != end && this->field(*tit, this->ghost_type).getSize() == 0; + ++tit) { + } type = *tit; - const array_type & vect = this->field(type, this->ghost_type); - UInt nb_data = this->getNbDataPerElem(type); - UInt nb_component = vect.getNbComponent(); - UInt size = (vect.getSize() * nb_component) / nb_data; - array_iterator it = vect.end_reinterpret(nb_data, size); + if (tit == end) + return this->end(); + + /// getting information for the field of the given type + const array_type & vect = this->field(type, this->ghost_type); + UInt nb_data_per_elem = this->getNbDataPerElem(type); + UInt nb_component = vect.getNbComponent(); + UInt size = (vect.getSize() * nb_component) / nb_data_per_elem; + + /// define element-wise iterator + array_iterator it = vect.begin_reinterpret(nb_data_per_elem, size); + array_iterator it_end = vect.end_reinterpret(nb_data_per_elem, size); + /// define data iterator + iterator rit = iterator(this->field, tit, end, it, it_end, + this->ghost_type, sub_element); + rit.setNbDataPerElem(this->nb_data_per_elem); + return rit; + } - iterator rit = - iterator(this->field, end, end, it, it, this->ghost_type, sub_element); - rit.setNbDataPerElem(this->nb_data_per_elem); - return rit; - } + virtual iterator end() { + field_type_iterator tit; + field_type_iterator end; + UInt sub_element = 0; + + tit = this->field.firstType(this->spatial_dimension, this->ghost_type, + this->element_kind); + end = this->field.lastType(this->spatial_dimension, this->ghost_type, + this->element_kind); + + ElementType type = *tit; + for (; tit != end; ++tit) + type = *tit; + + const array_type & vect = this->field(type, this->ghost_type); + UInt nb_data = this->getNbDataPerElem(type); + UInt nb_component = vect.getNbComponent(); + UInt size = (vect.getSize() * nb_component) / nb_data; + array_iterator it = vect.end_reinterpret(nb_data, size); + + iterator rit = iterator(this->field, end, end, it, it, this->ghost_type, + sub_element); + rit.setNbDataPerElem(this->nb_data_per_elem); + return rit; + } - /* ------------------------------------------------------------------------ */ - /* Class Members */ - /* ------------------------------------------------------------------------ */ + /* ------------------------------------------------------------------------ + */ + /* Class Members */ + /* ------------------------------------------------------------------------ + */ -protected: -}; + protected: + }; } // namespace dumpers } // namespace akantu #endif /* AKANTU_DUMPER_IGFEM_GENERIC_ELEMENTAL_FIELD_HH_ */ diff --git a/extra_packages/igfem/src/dumper_igfem_material_internal_field.hh b/extra_packages/igfem/src/dumper_igfem_material_internal_field.hh index 80c602966..5af3c8299 100644 --- a/extra_packages/igfem/src/dumper_igfem_material_internal_field.hh +++ b/extra_packages/igfem/src/dumper_igfem_material_internal_field.hh @@ -1,53 +1,58 @@ /** * @file dumper_igfem_material_internal_field.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief description of IGFEM material internal field * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ #ifndef AKANTU_DUMPER_IGFEM_MATERIAL_INTERNAL_FIELD_HH_ #define AKANTU_DUMPER_IGFEM_MATERIAL_INTERNAL_FIELD_HH_ /* -------------------------------------------------------------------------- */ #include "dumper_igfem_quadrature_points_field.hh" /* -------------------------------------------------------------------------- */ namespace akantu { namespace dumpers { -/* -------------------------------------------------------------------------- */ - -template -class IGFEMInternalMaterialField - : public IGFEMGenericElementalField, - igfem_quadrature_point_iterator> { - - /* ------------------------------------------------------------------------ */ - /* Typedefs */ - /* ------------------------------------------------------------------------ */ - -public: - typedef SingleType types; - typedef IGFEMGenericElementalField - parent; - typedef typename types::field_type field_type; - - /* ------------------------------------------------------------------------ */ - /* Constructors/Destructors */ - /* ------------------------------------------------------------------------ */ - - IGFEMInternalMaterialField(const field_type & field, - UInt spatial_dimension = _all_dimensions, - GhostType ghost_type = _not_ghost, - ElementKind kind = _ek_igfem) - : parent(field, spatial_dimension, ghost_type, kind) {} -}; + /* -------------------------------------------------------------------------- + */ + + template + class IGFEMInternalMaterialField + : public IGFEMGenericElementalField, + igfem_quadrature_point_iterator> { + + /* ------------------------------------------------------------------------ + */ + /* Typedefs */ + /* ------------------------------------------------------------------------ + */ + + public: + typedef SingleType types; + typedef IGFEMGenericElementalField + parent; + typedef typename types::field_type field_type; + + /* ------------------------------------------------------------------------ + */ + /* Constructors/Destructors */ + /* ------------------------------------------------------------------------ + */ + + IGFEMInternalMaterialField(const field_type & field, + UInt spatial_dimension = _all_dimensions, + GhostType ghost_type = _not_ghost, + ElementKind kind = _ek_igfem) + : parent(field, spatial_dimension, ghost_type, kind) {} + }; } // namespace dumpers } // namespace akantu #endif /* AKANTU_DUMPER_IGFEM_MATERIAL_INTERNAL_FIELD_HH_ */ diff --git a/extra_packages/igfem/src/dumper_igfem_quadrature_points_field.hh b/extra_packages/igfem/src/dumper_igfem_quadrature_points_field.hh index 38b060bba..ca9d90dce 100644 --- a/extra_packages/igfem/src/dumper_igfem_quadrature_points_field.hh +++ b/extra_packages/igfem/src/dumper_igfem_quadrature_points_field.hh @@ -1,140 +1,156 @@ /** * @file dumper_igfem_quadrature_points_field.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief description of IGFEM quadrature points field * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ #ifndef AKANTU_DUMPER_IGFEM_QUADRATURE_POINTS_FIELD_HH_ #define AKANTU_DUMPER_IGFEM_QUADRATURE_POINTS_FIELD_HH_ /* -------------------------------------------------------------------------- */ #include "dumper_igfem_elemental_field.hh" namespace akantu { namespace dumpers { -/* -------------------------------------------------------------------------- */ -template -class igfem_quadrature_point_iterator - : public igfem_element_iterator { - /* ------------------------------------------------------------------------ */ - /* Typedefs */ - /* ------------------------------------------------------------------------ */ -public: - typedef igfem_element_iterator - parent; - typedef typename types::data_type data_type; - typedef typename types::return_type return_type; - typedef typename types::field_type field_type; - typedef typename types::array_iterator array_iterator; - - /* ------------------------------------------------------------------------ */ - /* Constructors/Destructors */ - /* ------------------------------------------------------------------------ */ -public: - igfem_quadrature_point_iterator( - const field_type & field, const typename field_type::type_iterator & t_it, - const typename field_type::type_iterator & t_it_end, - const array_iterator & array_it, const array_iterator & array_it_end, - const GhostType ghost_type = _not_ghost, UInt sub_element = 0) - : parent(field, t_it, t_it_end, array_it, array_it_end, ghost_type, - sub_element) {} - - return_type operator*() { - const Vector & mat_internal_field = *this->array_it; - /// get nb data per sub element - UInt nb_sub_1_internal_points = - IGFEMHelper::getNbQuadraturePoints(*this->tit, 0); - UInt nb_sub_2_internal_points = - IGFEMHelper::getNbQuadraturePoints(*this->tit, 1); - - UInt nb_data = this->getNbDataPerElem(*(this->tit)) / - (nb_sub_1_internal_points + nb_sub_2_internal_points); - - UInt nb_sub_components = 0; - if (!(this->sub_element)) - nb_sub_components = nb_data * nb_sub_1_internal_points; - else - nb_sub_components = nb_data * nb_sub_2_internal_points; - - Vector sub_mat_internal_field(nb_sub_components); - - if (!(this->sub_element)) { - for (UInt i = 0; i < nb_sub_components; ++i) - sub_mat_internal_field(i) = mat_internal_field(i); - } else { - for (UInt i = 0; i < nb_sub_components; ++i) - sub_mat_internal_field(i) = - mat_internal_field(nb_data * nb_sub_1_internal_points + i); + /* -------------------------------------------------------------------------- + */ + template + class igfem_quadrature_point_iterator + : public igfem_element_iterator { + /* ------------------------------------------------------------------------ + */ + /* Typedefs */ + /* ------------------------------------------------------------------------ + */ + public: + typedef igfem_element_iterator + parent; + typedef typename types::data_type data_type; + typedef typename types::return_type return_type; + typedef typename types::field_type field_type; + typedef typename types::array_iterator array_iterator; + + /* ------------------------------------------------------------------------ + */ + /* Constructors/Destructors */ + /* ------------------------------------------------------------------------ + */ + public: + igfem_quadrature_point_iterator( + const field_type & field, + const typename field_type::type_iterator & t_it, + const typename field_type::type_iterator & t_it_end, + const array_iterator & array_it, const array_iterator & array_it_end, + const GhostType ghost_type = _not_ghost, UInt sub_element = 0) + : parent(field, t_it, t_it_end, array_it, array_it_end, ghost_type, + sub_element) {} + + return_type operator*() { + const Vector & mat_internal_field = *this->array_it; + /// get nb data per sub element + UInt nb_sub_1_internal_points = + IGFEMHelper::getNbQuadraturePoints(*this->tit, 0); + UInt nb_sub_2_internal_points = + IGFEMHelper::getNbQuadraturePoints(*this->tit, 1); + + UInt nb_data = this->getNbDataPerElem(*(this->tit)) / + (nb_sub_1_internal_points + nb_sub_2_internal_points); + + UInt nb_sub_components = 0; + if (!(this->sub_element)) + nb_sub_components = nb_data * nb_sub_1_internal_points; + else + nb_sub_components = nb_data * nb_sub_2_internal_points; + + Vector sub_mat_internal_field(nb_sub_components); + + if (!(this->sub_element)) { + for (UInt i = 0; i < nb_sub_components; ++i) + sub_mat_internal_field(i) = mat_internal_field(i); + } else { + for (UInt i = 0; i < nb_sub_components; ++i) + sub_mat_internal_field(i) = + mat_internal_field(nb_data * nb_sub_1_internal_points + i); + } + + return sub_mat_internal_field; } - - return sub_mat_internal_field; - } -}; - -// /* -------------------------------------------------------------------------- -// */ -// /* Fields type description */ -// /* -------------------------------------------------------------------------- -// */ -// template class iterator_type> -// class GenericQuadraturePointsField : -// public GenericElementalField { - -// public: - -// /* ------------------------------------------------------------------------ -// */ -// /* Typedefs */ -// /* ------------------------------------------------------------------------ -// */ - -// typedef iterator_type iterator; -// typedef typename types::field_type field_type; -// typedef typename iterator::it_type T; -// typedef GenericElementalField parent; - -// /* ------------------------------------------------------------------------ -// */ -// /* Constructors/Destructors */ -// /* ------------------------------------------------------------------------ -// */ - -// GenericQuadraturePointsField(const field_type & field, -// UInt spatial_dimension = _all_dimensions, -// GhostType ghost_type = _not_ghost, -// ElementKind element_kind = _ek_not_defined) : -// parent(field, spatial_dimension, ghost_type, element_kind) { } - -// /* ------------------------------------------------------------------------ -// */ -// /* Methods */ -// /* ------------------------------------------------------------------------ -// */ - -// virtual iterator begin() { -// iterator it = parent::begin(); -// return it; -// } - -// virtual iterator end () { -// iterator it = parent::end(); -// return it; -// } - -// }; - -/* -------------------------------------------------------------------------- */ + }; + + // /* + // -------------------------------------------------------------------------- + // */ + // /* Fields type description */ + // /* + // -------------------------------------------------------------------------- + // */ + // template class iterator_type> + // class GenericQuadraturePointsField : + // public GenericElementalField { + + // public: + + // /* + // ------------------------------------------------------------------------ + // */ + // /* Typedefs */ + // /* + // ------------------------------------------------------------------------ + // */ + + // typedef iterator_type iterator; + // typedef typename types::field_type field_type; + // typedef typename iterator::it_type T; + // typedef GenericElementalField parent; + + // /* + // ------------------------------------------------------------------------ + // */ + // /* Constructors/Destructors */ + // /* + // ------------------------------------------------------------------------ + // */ + + // GenericQuadraturePointsField(const field_type & field, + // UInt spatial_dimension = _all_dimensions, + // GhostType ghost_type = _not_ghost, + // ElementKind element_kind = _ek_not_defined) : + // parent(field, spatial_dimension, ghost_type, element_kind) { } + + // /* + // ------------------------------------------------------------------------ + // */ + // /* Methods */ + // /* + // ------------------------------------------------------------------------ + // */ + + // virtual iterator begin() { + // iterator it = parent::begin(); + // return it; + // } + + // virtual iterator end () { + // iterator it = parent::end(); + // return it; + // } + + // }; + + /* -------------------------------------------------------------------------- + */ } // namespace dumpers } // namespace akantu #endif /* AKANTU_DUMPER_IGFEM_QUADRATURE_POINTS_FIELD_HH_ */ diff --git a/extra_packages/igfem/src/element_class_igfem.hh b/extra_packages/igfem/src/element_class_igfem.hh index 69b7ca250..c98cc6cf6 100644 --- a/extra_packages/igfem/src/element_class_igfem.hh +++ b/extra_packages/igfem/src/element_class_igfem.hh @@ -1,297 +1,304 @@ /** * @file element_class_igfem.hh * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * * @brief Specialization for interface-enriched finite elements * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_ELEMENT_CLASS_IGFEM_HH_ #define AKANTU_ELEMENT_CLASS_IGFEM_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ template class InterpolationElement { public: using interpolation_property = InterpolationProperty; /* ------------------------------------------------------------------------ */ /* Member functions */ /* ------------------------------------------------------------------------ */ public: static void assembleShapes(const Vector & parent_interpolation, const Vector & sub_interpolation, Vector & interpolation, UInt sub_element = 0) { /// N1, N2, N3 of parent triangle UInt nb_nodes_parent = InterpolationElement< interpolation_property::parent_interpolation_type>::getShapeSize(); for (UInt i = 0; i < nb_nodes_parent; ++i) { interpolation(i) = parent_interpolation(i); } /// add the enrichment UInt * enriched_node = enrichments[sub_element]; for (UInt e = 0; e < nb_enrichments; ++e) { interpolation(nb_nodes_parent + e) = sub_interpolation(enriched_node[e]); } } static void assembleShapeDerivatives(const Matrix & parent_interpolation, const Matrix & sub_interpolation, Matrix & interpolation, UInt sub_element = 0) { /// N1, N2, N3 of parent triangle UInt nb_nodes_parent = InterpolationElement< interpolation_property::parent_interpolation_type>::getShapeSize(); for (UInt i = 0; i < nb_nodes_parent; ++i) { Vector ip(interpolation(i)); ip = parent_interpolation(i); } /// add the enrichment UInt * enriched_node = enrichments[sub_element]; for (UInt e = 0; e < nb_enrichments; ++e) { Vector ip(interpolation(nb_nodes_parent + e)); ip = sub_interpolation(enriched_node[e]); } } static void interpolate(const Matrix & nodal_values, const Vector & shapes, Vector & interpolated) { Matrix interpm(interpolated.storage(), nodal_values.rows(), 1); Matrix shapesm(shapes.storage(), interpolation_property::nb_nodes_per_element, 1); interpm.mul(nodal_values, shapesm); } public: static AKANTU_GET_MACRO_NOT_CONST( ShapeSize, interpolation_property::nb_nodes_per_element, UInt); static AKANTU_GET_MACRO_NOT_CONST( ShapeDerivativesSize, (interpolation_property::nb_nodes_per_element * interpolation_property::natural_space_dimension), UInt); static AKANTU_GET_MACRO_NOT_CONST( NaturalSpaceDimension, interpolation_property::natural_space_dimension, UInt); static AKANTU_GET_MACRO_NOT_CONST( NbNodesPerInterpolationElement, interpolation_property::nb_nodes_per_element, UInt); - static AKANTU_GET_MACRO_NOT_CONST(NbSubElements, interpolation_property::nb_sub_elements, UInt); + static AKANTU_GET_MACRO_NOT_CONST(NbSubElements, + interpolation_property::nb_sub_elements, + UInt); static UInt * getSubElementConnectivity(UInt t = 0) { return &(interpolation_property::sub_element_connectivity[t]); }; - static UInt getNbEnrichments() { return interpolation_property::nb_enrichments; }; - static UInt * getSubElementEnrichments(UInt t = 0) { return &(interpolation_property::enrichments[t]); }; + static UInt getNbEnrichments() { + return interpolation_property::nb_enrichments; + }; + static UInt * getSubElementEnrichments(UInt t = 0) { + return &(interpolation_property::enrichments[t]); + }; protected: /// storage of the subelement local connectivity static UInt sub_element_connectivity_vect[]; /// local connectivity of subelements static UInt * sub_element_connectivity[]; /// nb of subelements static UInt nb_sub_elements; /// storage of enrichments static UInt enrichment_vect[]; static UInt * enrichments[]; static UInt nb_enrichments; }; } // namespace akantu #include "interpolation_element_igfem_tmpl.hh" namespace akantu { /* -------------------------------------------------------------------------- */ #define AKANTU_DEFINE_IGFEM_ELEMENT_CLASS_PROPERTY( \ elem_type, geom_type, interp_type, parent_el_type, sub_el_type_1, \ sub_el_type_2, elem_kind, sp, min_int_order) \ template <> struct ElementClassProperty { \ static const GeometricalType geometrical_type{geom_type}; \ static const InterpolationType interpolation_type{interp_type}; \ static const ElementType parent_element_type{parent_el_type}; \ static const ElementType sub_element_type_1{sub_el_type_1}; \ static const ElementType sub_element_type_2{sub_el_type_2}; \ static const ElementKind element_kind{elem_kind}; \ static const UInt spatial_dimension{sp}; \ static const UInt minimal_integration_order{min_int_order}; \ } /* -------------------------------------------------------------------------- */ template class ElementClass : public GeometricalElement< ElementClassProperty::geometrical_type>, public InterpolationElement< ElementClassProperty::interpolation_type> { protected: using geometrical_element = GeometricalElement::geometrical_type>; using interpolation_element = InterpolationElement< ElementClassProperty::interpolation_type>; using parent_element = ElementClass::parent_element_type>; using element_property = ElementClassProperty; using interpolation_property = typename interpolation_element::interpolation_property; /* ------------------------------------------------------------------------ */ /* Member functions */ /* ------------------------------------------------------------------------ */ public: static void getSubElementCoords(const Matrix & element_coords, Matrix & sub_coords, const UInt sub_element) { /// get the sub_element_type /// constexrp ElementType sub_el_type = getSubElementType(sub_element); UInt nb_nodes_sub_el = 0; switch (sub_element) { case 0: nb_nodes_sub_el = ElementClass::sub_element_type_1>:: getNbNodesPerInterpolationElement(); break; case 1: nb_nodes_sub_el = ElementClass::sub_element_type_2>:: getNbNodesPerInterpolationElement(); break; } for (UInt i = 0; i < nb_nodes_sub_el; ++i) { UInt lc = InterpolationElement< ElementClassProperty::interpolation_type>:: sub_element_connectivity[sub_element][i]; Vector sub_c(sub_coords(i)); sub_c = element_coords(lc); } } static void getParentCoords(const Matrix & element_coords, Matrix & parent_coords) { const ElementType parent_type = ElementClassProperty::parent_element_type; UInt nb_nodes_parent_el = ElementClass::getNbNodesPerInterpolationElement(); for (UInt i = 0; i < nb_nodes_parent_el; ++i) { Vector parent_c(parent_coords(i)); parent_c = element_coords(i); } } /// map the points from the reference domain of the subelement to the physical /// domain static void mapToPhysicalDomain(const Matrix & element_coords, Matrix & sub_coords, Matrix & sub_shapes, Matrix & physical_points, UInt sub_element = 0) { /// get the sub_element_type getSubElementCoords(element_coords, sub_coords, sub_element); /// map the points of the subelements in the physical domain switch (sub_element) { case 0: ElementClass::sub_element_type_1>:: interpolate(sub_coords, sub_shapes, physical_points); break; case 1: ElementClass::sub_element_type_2>:: interpolate(sub_coords, sub_shapes, physical_points); break; } } /// map the points from the physical domain to the parent reference domain static void mapToParentRefDomain(const Matrix & element_coords, Matrix & parent_coords, Matrix & physical_points, Matrix & natural_coords) { const ElementType parent_type = ElementClassProperty::parent_element_type; getParentCoords(element_coords, parent_coords); /// map the points from the physical domain into the parent reference domain ElementClass::inverseMap(physical_points, parent_coords, natural_coords); } /// map the points from the subelement reference domain to the parent /// reference domain static void mapFromSubRefToParentRef(const Matrix & element_coords, Matrix & parent_coords, Matrix & sub_coords, Matrix & sub_shapes, Matrix & physical_points, Matrix & natural_points, UInt /*nb_points*/, UInt sub_element) { mapToPhysicalDomain(element_coords, sub_coords, sub_shapes, physical_points, sub_element); mapToParentRefDomain(element_coords, parent_coords, physical_points, natural_points); } static void mapFromSubRefToParentRef(const Matrix & element_coords, Matrix & sub_coords, Matrix & parent_coords, Matrix & sub_shapes, Matrix & physical_points, Matrix & parent_el_natural_coords, UInt sub_element) { mapToPhysicalDomain(element_coords, sub_coords, sub_shapes, physical_points, sub_element); mapToParentRefDomain(element_coords, parent_coords, physical_points, parent_el_natural_coords); } /// compute the normal of a surface defined by the function f static inline void computeNormalsOnNaturalCoordinates(const Matrix & /*coord*/, Matrix & /*f*/, Matrix & /*normals*/) { AKANTU_TO_IMPLEMENT(); } /// determine orientation of the element with respect to the interface static inline UInt getOrientation(const Vector & is_inside); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: static AKANTU_GET_MACRO_NOT_CONST(Kind, _ek_igfem, ElementKind); static ElementType getP1ElementType() { AKANTU_TO_IMPLEMENT(); }; static AKANTU_GET_MACRO_NOT_CONST( SpatialDimension, ElementClassProperty::spatial_dimension, UInt); static ElementType & getFacetType(UInt /*t*/ = 0) { AKANTU_TO_IMPLEMENT(); } static ElementType * getFacetTypeInternal() { AKANTU_TO_IMPLEMENT(); } + private: }; } // namespace akantu /* -------------------------------------------------------------------------- */ #include "geometrical_element_igfem.hh" /* -------------------------------------------------------------------------- */ #include "element_class_igfem_segment_3_inline_impl.hh" #include "element_class_igfem_triangle_4_inline_impl.hh" #include "element_class_igfem_triangle_5_inline_impl.hh" /* -------------------------------------------------------------------------- */ #endif /* AKANTU_ELEMENT_CLASS_IGFEM_HH_ */ diff --git a/extra_packages/igfem/src/fe_engine_template_tmpl_igfem.hh b/extra_packages/igfem/src/fe_engine_template_tmpl_igfem.hh index 3fe2d31a8..276ab296d 100644 --- a/extra_packages/igfem/src/fe_engine_template_tmpl_igfem.hh +++ b/extra_packages/igfem/src/fe_engine_template_tmpl_igfem.hh @@ -1,117 +1,115 @@ /** * @file shape_igfem_inline_impl.hh * * @author Aurelia Isabel Cuba Ramos * * @brief ShapeIGFEM inline implementation * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "integrator_gauss_igfem.hh" #include "shape_igfem.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_FE_ENGINE_TEMPLATE_TMPL_IGFEM_HH_ #define AKANTU_FE_ENGINE_TEMPLATE_TMPL_IGFEM_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ /* compatibility functions */ /* -------------------------------------------------------------------------- */ template <> inline void FEEngineTemplate:: - initShapeFunctions(const Array & nodes, - GhostType ghost_type) { + initShapeFunctions(const Array & nodes, GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh::type_iterator it = mesh.firstType(element_dimension, ghost_type, _ek_igfem); Mesh::type_iterator end = mesh.lastType(element_dimension, ghost_type, _ek_igfem); for (; it != end; ++it) { ElementType type = *it; integrator.initIntegrator(nodes, type, ghost_type); #define INIT(_type) \ do { \ const Matrix & all_quads = \ integrator.getIntegrationPoints<_type>(ghost_type); \ const Matrix & quads_1 = integrator.getIntegrationPoints< \ ElementClassProperty<_type>::sub_element_type_1>(ghost_type); \ const Matrix & quads_2 = integrator.getIntegrationPoints< \ ElementClassProperty<_type>::sub_element_type_2>(ghost_type); \ shape_functions.initShapeFunctions(nodes, all_quads, quads_1, quads_2, \ _type, ghost_type); \ } while (0) AKANTU_BOOST_IGFEM_ELEMENT_SWITCH(INIT); #undef INIT } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template <> inline void FEEngineTemplate:: computeIntegrationPointsCoordinates( Array & quadrature_points_coordinates, ElementType type, - GhostType ghost_type, - const Array & filter_elements) const { + GhostType ghost_type, const Array & filter_elements) const { const Array & nodes_coordinates = mesh.getNodes(); UInt spatial_dimension = mesh.getSpatialDimension(); /// create an array with the nodal coordinates that need to be /// interpolated. The nodal coordinates of the enriched nodes need /// to be set to zero, because they represent the enrichment of the /// position field, and the enrichments for this field are all zero! /// There is no gap in the mesh! Array igfem_nodes(nodes_coordinates.getSize(), spatial_dimension); shape_functions.extractValuesAtStandardNodes(nodes_coordinates, igfem_nodes, ghost_type); interpolateOnIntegrationPoints(igfem_nodes, quadrature_points_coordinates, spatial_dimension, type, ghost_type, filter_elements); } /* -------------------------------------------------------------------------- */ template <> inline void FEEngineTemplate:: computeIntegrationPointsCoordinates( ElementTypeMapArray & quadrature_points_coordinates, const ElementTypeMapArray * filter_elements) const { const Array & nodes_coordinates = mesh.getNodes(); UInt spatial_dimension = mesh.getSpatialDimension(); /// create an array with the nodal coordinates that need to be /// interpolated. The nodal coordinates of the enriched nodes need /// to be set to zero, because they represent the enrichment of the /// position field, and the enrichments for this field are all zero! /// There is no gap in the mesh! Array igfem_nodes(nodes_coordinates.getSize(), spatial_dimension); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; shape_functions.extractValuesAtStandardNodes(nodes_coordinates, igfem_nodes, ghost_type); } interpolateOnIntegrationPoints(igfem_nodes, quadrature_points_coordinates, filter_elements); } } // namespace akantu #endif /* AKANTU_FE_ENGINE_TEMPLATE_TMPL_IGFEM_HH_ */ diff --git a/extra_packages/igfem/src/igfem_helper.hh b/extra_packages/igfem/src/igfem_helper.hh index 687e9653b..01a103fa5 100644 --- a/extra_packages/igfem/src/igfem_helper.hh +++ b/extra_packages/igfem/src/igfem_helper.hh @@ -1,148 +1,147 @@ /** * @file dumper_igfem_element_iterator.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief Helper class to return sub element information * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ #ifndef AKANTU_IGFEM_HELPER_HH_ #define AKANTU_IGFEM_HELPER_HH_ /* -------------------------------------------------------------------------- */ #include "element_class.hh" /* -------------------------------------------------------------------------- */ namespace akantu { class FEEngine; template struct ElementTypeIGFEMData { static ElementType igfem_element_types[]; static UInt nb_igfem_types; }; struct IGFEMHelper { template static VectorProxy getIGFEMElementTypes() { return VectorProxy( ElementTypeIGFEMData::igfem_element_types, ElementTypeIGFEMData::nb_igfem_types); } /// get the number of nodes for a given sub-element static UInt getNbNodesPerSubElement(ElementType type, const UInt sub_element) { UInt nb_nodes_per_sub_element = 0; #define GET_NB_NODES_PER_SUB_ELEMENT(type) \ switch (sub_element) { \ case 0: \ nb_nodes_per_sub_element = \ ElementClass::sub_element_type_1>:: \ getNbNodesPerInterpolationElement(); \ break; \ case 1: \ nb_nodes_per_sub_element = \ ElementClass::sub_element_type_2>:: \ getNbNodesPerInterpolationElement(); \ break; \ } AKANTU_BOOST_IGFEM_ELEMENT_SWITCH(GET_NB_NODES_PER_SUB_ELEMENT); #undef GET_NB_NODES_PER_SUB_ELEMENT return nb_nodes_per_sub_element; } /// get the connectivity for a given sub-element static UInt * getSubElementConnectivity(ElementType type, const UInt sub_element) { UInt * sub_element_connectivity = NULL; #define GET_SUB_ELEMENT_CONNECTIVITY(type) \ sub_element_connectivity = \ ElementClass::getSubElementConnectivity(sub_element); AKANTU_BOOST_IGFEM_ELEMENT_SWITCH(GET_SUB_ELEMENT_CONNECTIVITY); #undef GET_SUB_ELEMENT_CONNECTIVITY return sub_element_connectivity; } /// get the sub-element type static ElementType getSubElementType(ElementType type, const UInt sub_element) { ElementType sub_type = _not_defined; #define GET_SUB_ELEMENT_TYPE(type) \ switch (sub_element) { \ case 0: \ sub_type = ElementClassProperty::sub_element_type_1; \ break; \ case 1: \ sub_type = ElementClassProperty::sub_element_type_2; \ break; \ } AKANTU_BOOST_IGFEM_ELEMENT_SWITCH(GET_SUB_ELEMENT_TYPE); #undef GET_SUB_ELEMENT_TYPE return sub_type; } /// get the nb of quads for one sub element type - static UInt getNbQuadraturePoints(ElementType type, - const UInt sub_element) { + static UInt getNbQuadraturePoints(ElementType type, const UInt sub_element) { UInt nb_quad_points = 0; #define GET_NB_QUADS(type) \ switch (sub_element) { \ case 0: \ nb_quad_points = GaussIntegrationElement::sub_element_type_1>::getNbQuadraturePoints(); \ break; \ case 1: \ nb_quad_points = GaussIntegrationElement::sub_element_type_2>::getNbQuadraturePoints(); \ break; \ } AKANTU_BOOST_IGFEM_ELEMENT_SWITCH(GET_NB_QUADS); #undef GET_NB_QUADS return nb_quad_points; } /// get the nb of parent nodes of a given igfem element type static UInt getNbParentNodes(ElementType type) { UInt nb_parent_nodes = 0; #define GET_NB_PARENT_NODES(type) \ nb_parent_nodes = \ ElementClass::parent_element_type>:: \ getNbNodesPerInterpolationElement(); AKANTU_BOOST_IGFEM_ELEMENT_SWITCH(GET_NB_PARENT_NODES); #undef GET_NB_PARENT_NODES return nb_parent_nodes; } /// get the nb of parent nodes of a given igfem element type static UInt getNbEnrichedNodes(ElementType type) { UInt nb_enriched_nodes = 0; #define GET_NB_ENRICHED_NODES(type) \ nb_enriched_nodes = ElementClass::getNbEnrichments(); AKANTU_BOOST_IGFEM_ELEMENT_SWITCH(GET_NB_ENRICHED_NODES); #undef GET_NB_ENRICHED_NODES return nb_enriched_nodes; } /// get the nb of quads for one sub element type static UInt getElementOrientation(ElementType type, const Vector & is_inside) { UInt orientation = 0; #define GET_ORIENTATION(type) \ orientation = ElementClass::getOrientation(is_inside); AKANTU_BOOST_IGFEM_ELEMENT_SWITCH(GET_ORIENTATION); #undef GET_ORIENTATION return orientation; } }; } // namespace akantu #endif /* AKANTU_IGFEM_HELPER_HH_ */ diff --git a/extra_packages/igfem/src/integrator_gauss_igfem.hh b/extra_packages/igfem/src/integrator_gauss_igfem.hh index c32092989..02f75885d 100644 --- a/extra_packages/igfem/src/integrator_gauss_igfem.hh +++ b/extra_packages/igfem/src/integrator_gauss_igfem.hh @@ -1,123 +1,119 @@ /** * @file integrator_gauss_igfem.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief Gauss integration facilities for IGFEM * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #ifndef AKANTU_INTEGRATOR_IGFEM_HH_ #define AKANTU_INTEGRATOR_IGFEM_HH_ /* -------------------------------------------------------------------------- */ #include "integrator.hh" /* -------------------------------------------------------------------------- */ namespace akantu { template class IntegratorGauss<_ek_igfem, IOF> : public Integrator { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: IntegratorGauss(const Mesh & mesh, const ID & id = "integrator_gauss"); virtual ~IntegratorGauss(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: - inline void initIntegrator(const Array & nodes, - ElementType type, + inline void initIntegrator(const Array & nodes, ElementType type, GhostType ghost_type); /// precompute jacobians on elements of type "type" template void precomputeJacobiansOnQuadraturePoints(const Array & nodes, GhostType ghost_type); /// integrate f on the element "elem" of type "type" template inline void integrateOnElement(const Array & f, Real * intf, UInt nb_degree_of_freedom, UInt elem, GhostType ghost_type) const; /// integrate f for all elements of type "type" template void integrate(const Array & in_f, Array & intf, UInt nb_degree_of_freedom, GhostType ghost_type, const Array & filter_elements) const; /// integrate one element scalar value on all elements of type "type" template Real integrate(const Vector & in_f, UInt index, GhostType ghost_type) const; /// integrate scalar field in_f template Real integrate(const Array & in_f, GhostType ghost_type, const Array & filter_elements) const; /// integrate partially around a quadrature point (@f$ intf_q = f_q * J_q * /// w_q @f$) template - void integrateOnIntegrationPoints(const Array & in_f, - Array & intf, - UInt nb_degree_of_freedom, - GhostType ghost_type, - const Array & filter_elements) const; + void + integrateOnIntegrationPoints(const Array & in_f, Array & intf, + UInt nb_degree_of_freedom, GhostType ghost_type, + const Array & filter_elements) const; /// return a vector with quadrature points natural coordinates template const Matrix & getIntegrationPoints(GhostType ghost_type) const; /// return the number of quadrature points for a given element type template - inline UInt - getNbIntegrationPoints(GhostType ghost_type = _not_ghost) const; + inline UInt getNbIntegrationPoints(GhostType ghost_type = _not_ghost) const; /// compute the vector of quadrature points natural coordinates template void computeQuadraturePoints(GhostType ghost_type); /// check that the jacobians are not negative - template - void checkJacobians(GhostType ghost_type) const; + template void checkJacobians(GhostType ghost_type) const; public: /// compute the jacobians on quad points for a given element template void computeJacobianOnQuadPointsByElement(const Matrix & node_coords, Vector & jacobians); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: inline void integrate(Real * f, Real * jac, Real * inte, UInt nb_degree_of_freedom, UInt nb_quadrature_points) const; private: ElementTypeMap> quadrature_points; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "integrator_gauss_igfem_inline_impl.hh" } // namespace akantu #endif /*AKANTU_INTEGRATOR_IGFEM_HH_ */ /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ diff --git a/extra_packages/igfem/src/integrator_gauss_igfem_inline_impl.hh b/extra_packages/igfem/src/integrator_gauss_igfem_inline_impl.hh index 376ccba79..6d182f7e6 100644 --- a/extra_packages/igfem/src/integrator_gauss_igfem_inline_impl.hh +++ b/extra_packages/igfem/src/integrator_gauss_igfem_inline_impl.hh @@ -1,451 +1,449 @@ /** * @file integrator_gauss_igfem.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief Inline functions of gauss integrator for the case of IGFEM * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ } // namespace akantu #include "fe_engine.hh" #if defined(AKANTU_DEBUG_TOOLS) #include "aka_debug_tools.hh" #endif namespace akantu { /* -------------------------------------------------------------------------- */ #define INIT_INTEGRATOR(type) \ computeQuadraturePoints(ghost_type); \ precomputeJacobiansOnQuadraturePoints(nodes, ghost_type); \ checkJacobians(ghost_type); template -inline void -IntegratorGauss<_ek_igfem, IOF>::initIntegrator(const Array & nodes, - ElementType type, - GhostType ghost_type) { +inline void IntegratorGauss<_ek_igfem, IOF>::initIntegrator( + const Array & nodes, ElementType type, GhostType ghost_type) { AKANTU_BOOST_IGFEM_ELEMENT_SWITCH(INIT_INTEGRATOR); } #undef INIT_INTEGRATOR /* -------------------------------------------------------------------------- */ template template -inline UInt IntegratorGauss<_ek_igfem, IOF>::getNbIntegrationPoints( - GhostType) const { +inline UInt +IntegratorGauss<_ek_igfem, IOF>::getNbIntegrationPoints(GhostType) const { const ElementType sub_type_1 = ElementClassProperty::sub_element_type_1; const ElementType sub_type_2 = ElementClassProperty::sub_element_type_2; UInt nb_quad_points_sub_1 = GaussIntegrationElement::getNbQuadraturePoints(); UInt nb_quad_points_sub_2 = GaussIntegrationElement::getNbQuadraturePoints(); return (nb_quad_points_sub_1 + nb_quad_points_sub_2); } /* -------------------------------------------------------------------------- */ template template inline void IntegratorGauss<_ek_igfem, IOF>::integrateOnElement( const Array & f, Real * intf, UInt nb_degree_of_freedom, const UInt elem, GhostType ghost_type) const { Array & jac_loc = jacobians(type, ghost_type); UInt nb_quadrature_points = getNbIntegrationPoints(); AKANTU_DEBUG_ASSERT(f.getNbComponent() == nb_degree_of_freedom, "The vector f do not have the good number of component."); Real * f_val = f.storage() + elem * f.getNbComponent(); Real * jac_val = jac_loc.storage() + elem * nb_quadrature_points; integrate(f_val, jac_val, intf, nb_degree_of_freedom, nb_quadrature_points); } /* -------------------------------------------------------------------------- */ template template inline Real IntegratorGauss<_ek_igfem, IOF>::integrate( const Vector & in_f, UInt index, GhostType ghost_type) const { const Array & jac_loc = jacobians(type, ghost_type); UInt nb_quadrature_points = getNbIntegrationPoints(); AKANTU_DEBUG_ASSERT(in_f.size() == nb_quadrature_points, "The vector f do not have nb_quadrature_points entries."); Real * jac_val = jac_loc.storage() + index * nb_quadrature_points; Real intf; integrate(in_f.storage(), jac_val, &intf, 1, nb_quadrature_points); return intf; return 0.; } /* -------------------------------------------------------------------------- */ template inline void IntegratorGauss<_ek_igfem, IOF>::integrate(Real * f, Real * jac, Real * inte, UInt nb_degree_of_freedom, UInt nb_quadrature_points) const { memset(inte, 0, nb_degree_of_freedom * sizeof(Real)); Real * cjac = jac; for (UInt q = 0; q < nb_quadrature_points; ++q) { for (UInt dof = 0; dof < nb_degree_of_freedom; ++dof) { inte[dof] += *f * *cjac; ++f; } ++cjac; } } /* -------------------------------------------------------------------------- */ template template inline const Matrix & IntegratorGauss<_ek_igfem, IOF>::getIntegrationPoints( GhostType ghost_type) const { AKANTU_DEBUG_ASSERT( quadrature_points.exists(type, ghost_type), "Quadrature points for type " << quadrature_points.printType(type, ghost_type) << " have not been initialized." << " Did you use 'computeQuadraturePoints' function ?"); return quadrature_points(type, ghost_type); } /* -------------------------------------------------------------------------- */ template template -inline void IntegratorGauss<_ek_igfem, IOF>::computeQuadraturePoints( - GhostType ghost_type) { +inline void +IntegratorGauss<_ek_igfem, IOF>::computeQuadraturePoints(GhostType ghost_type) { /// typedef for the two subelement_types and the parent element type const ElementType sub_type_1 = ElementClassProperty::sub_element_type_1; const ElementType sub_type_2 = ElementClassProperty::sub_element_type_2; /// store the quadrature points on the two subelements Matrix & quads_sub_1 = quadrature_points(sub_type_1, ghost_type); Matrix & quads_sub_2 = quadrature_points(sub_type_2, ghost_type); quads_sub_1 = GaussIntegrationElement::getQuadraturePoints(); quads_sub_2 = GaussIntegrationElement::getQuadraturePoints(); /// store all quad points for the current type UInt nb_quad_points_sub_1 = GaussIntegrationElement::getNbQuadraturePoints(); UInt nb_quad_points_sub_2 = GaussIntegrationElement::getNbQuadraturePoints(); UInt spatial_dimension = mesh.getSpatialDimension(); Matrix & quads = quadrature_points(type, ghost_type); quads = Matrix(spatial_dimension, nb_quad_points_sub_1 + nb_quad_points_sub_2); Matrix quads_1(quads.storage(), quads.rows(), nb_quad_points_sub_1); quads_1 = quads_sub_1; Matrix quads_2(quads.storage() + quads.rows() * nb_quad_points_sub_1, quads.rows(), nb_quad_points_sub_2); quads_2 = quads_sub_2; } /* -------------------------------------------------------------------------- */ template template inline void IntegratorGauss<_ek_igfem, IOF>::computeJacobianOnQuadPointsByElement( const Matrix & node_coords, Vector & jacobians) { /// optimize: get the matrix from the ElementTypeMap Matrix quad = GaussIntegrationElement::getQuadraturePoints(); // jacobian ElementClass::computeJacobian(quad, node_coords, jacobians); } /* -------------------------------------------------------------------------- */ template -inline IntegratorGauss<_ek_igfem, IOF>::IntegratorGauss( - const Mesh & mesh, const ID & id) +inline IntegratorGauss<_ek_igfem, IOF>::IntegratorGauss(const Mesh & mesh, + const ID & id) : Integrator(mesh, id) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template -inline void IntegratorGauss<_ek_igfem, IOF>::checkJacobians( - GhostType ghost_type) const { +inline void +IntegratorGauss<_ek_igfem, IOF>::checkJacobians(GhostType ghost_type) const { AKANTU_DEBUG_IN(); /// typedef for the two subelement_types and the parent element type const ElementType sub_type_1 = ElementClassProperty::sub_element_type_1; const ElementType sub_type_2 = ElementClassProperty::sub_element_type_2; UInt nb_quad_points_sub_1 = GaussIntegrationElement::getNbQuadraturePoints(); UInt nb_quad_points_sub_2 = GaussIntegrationElement::getNbQuadraturePoints(); UInt nb_quadrature_points = nb_quad_points_sub_1 + nb_quad_points_sub_2; UInt nb_element; nb_element = mesh.getConnectivity(type, ghost_type).getSize(); Real * jacobians_val = jacobians(type, ghost_type).storage(); for (UInt i = 0; i < nb_element * nb_quadrature_points; ++i, ++jacobians_val) { if (*jacobians_val < 0) AKANTU_ERROR( "Negative jacobian computed," << " possible problem in the element node ordering (Quadrature Point " << i % nb_quadrature_points << ":" << i / nb_quadrature_points << ":" << type << ":" << ghost_type << ")"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template inline void IntegratorGauss<_ek_igfem, IOF>::precomputeJacobiansOnQuadraturePoints( const Array & nodes, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// typedef for the two subelement_types and the parent element type const ElementType sub_type_1 = ElementClassProperty::sub_element_type_1; const ElementType sub_type_2 = ElementClassProperty::sub_element_type_2; UInt spatial_dimension = mesh.getSpatialDimension(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); /// get the number of nodes for the subelements and the parent element UInt nb_nodes_sub_1 = ElementClass::getNbNodesPerInterpolationElement(); UInt nb_nodes_sub_2 = ElementClass::getNbNodesPerInterpolationElement(); UInt nb_quadrature_points_sub_1 = GaussIntegrationElement::getNbQuadraturePoints(); UInt nb_quadrature_points_sub_2 = GaussIntegrationElement::getNbQuadraturePoints(); UInt nb_quadrature_points = nb_quadrature_points_sub_1 + nb_quadrature_points_sub_2; UInt nb_element = mesh.getNbElement(type, ghost_type); Array * jacobians_tmp; if (!jacobians.exists(type, ghost_type)) jacobians_tmp = &jacobians.alloc(nb_element * nb_quadrature_points, 1, type, ghost_type); else { jacobians_tmp = &jacobians(type, ghost_type); jacobians_tmp->resize(nb_element * nb_quadrature_points); } Array::vector_iterator jacobians_it = jacobians_tmp->begin_reinterpret(nb_quadrature_points, nb_element); Vector weights_sub_1 = GaussIntegrationElement::getWeights(); Vector weights_sub_2 = GaussIntegrationElement::getWeights(); Array x_el(0, spatial_dimension * nb_nodes_per_element); FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type); Array::const_matrix_iterator x_it = x_el.begin(spatial_dimension, nb_nodes_per_element); // Matrix local_coord(spatial_dimension, nb_nodes_per_element); for (UInt elem = 0; elem < nb_element; ++elem, ++jacobians_it, ++x_it) { const Matrix & X = *x_it; Matrix sub_1_coords(spatial_dimension, nb_nodes_sub_1); Matrix sub_2_coords(spatial_dimension, nb_nodes_sub_2); ElementClass::getSubElementCoords(X, sub_1_coords, 0); ElementClass::getSubElementCoords(X, sub_2_coords, 1); Vector & J = *jacobians_it; /// initialize vectors to store the jacobians for each subelement Vector J_sub_1(nb_quadrature_points_sub_1); Vector J_sub_2(nb_quadrature_points_sub_2); computeJacobianOnQuadPointsByElement(sub_1_coords, J_sub_1); computeJacobianOnQuadPointsByElement(sub_2_coords, J_sub_2); J_sub_1 *= weights_sub_1; J_sub_2 *= weights_sub_2; /// copy results into the jacobian vector for this element for (UInt i = 0; i < nb_quadrature_points_sub_1; ++i) { J(i) = J_sub_1(i); } for (UInt i = 0; i < nb_quadrature_points_sub_2; ++i) { J(i + nb_quadrature_points_sub_1) = J_sub_2(i); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template inline void IntegratorGauss<_ek_igfem, IOF>::integrate( const Array & in_f, Array & intf, UInt nb_degree_of_freedom, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(jacobians.exists(type, ghost_type), "No jacobians for the type " << jacobians.printType(type, ghost_type)); const Matrix & quads = quadrature_points(type, ghost_type); UInt nb_points = quads.cols(); const Array & jac_loc = jacobians(type, ghost_type); Array::const_matrix_iterator J_it; Array::matrix_iterator inte_it; Array::const_matrix_iterator f_it; UInt nb_element; Array * filtered_J = NULL; if (filter_elements != empty_filter) { nb_element = filter_elements.getSize(); filtered_J = new Array(0, jac_loc.getNbComponent()); FEEngine::filterElementalData(mesh, jac_loc, *filtered_J, type, ghost_type, filter_elements); const Array & cfiltered_J = *filtered_J; // \todo temporary patch J_it = cfiltered_J.begin_reinterpret(nb_points, 1, nb_element); } else { nb_element = mesh.getNbElement(type, ghost_type); J_it = jac_loc.begin_reinterpret(nb_points, 1, nb_element); } intf.resize(nb_element); f_it = in_f.begin_reinterpret(nb_degree_of_freedom, nb_points, nb_element); inte_it = intf.begin_reinterpret(nb_degree_of_freedom, 1, nb_element); for (UInt el = 0; el < nb_element; ++el, ++J_it, ++f_it, ++inte_it) { const Matrix & f = *f_it; const Matrix & J = *J_it; Matrix & inte_f = *inte_it; inte_f.mul(f, J); } delete filtered_J; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template inline Real IntegratorGauss<_ek_igfem, IOF>::integrate( const Array & in_f, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(jacobians.exists(type, ghost_type), "No jacobians for the type " << jacobians.printType(type, ghost_type)); Array intfv(0, 1); integrate(in_f, intfv, 1, ghost_type, filter_elements); UInt nb_values = intfv.getSize(); if (nb_values == 0) return 0.; UInt nb_values_to_sum = nb_values >> 1; std::sort(intfv.begin(), intfv.end()); // as long as the half is not empty while (nb_values_to_sum) { UInt remaining = (nb_values - 2 * nb_values_to_sum); if (remaining) intfv(nb_values - 2) += intfv(nb_values - 1); // sum to consecutive values and store the sum in the first half for (UInt i = 0; i < nb_values_to_sum; ++i) { intfv(i) = intfv(2 * i) + intfv(2 * i + 1); } nb_values = nb_values_to_sum; nb_values_to_sum >>= 1; } AKANTU_DEBUG_OUT(); return intfv(0); } /* -------------------------------------------------------------------------- */ template template inline void IntegratorGauss<_ek_igfem, IOF>::integrateOnIntegrationPoints( const Array & in_f, Array & intf, UInt nb_degree_of_freedom, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(jacobians.exists(type, ghost_type), "No jacobians for the type " << jacobians.printType(type, ghost_type)); UInt nb_element; const Matrix & quads = quadrature_points(type, ghost_type); UInt nb_points = quads.cols(); const Array & jac_loc = jacobians(type, ghost_type); Array::const_scalar_iterator J_it; Array::vector_iterator inte_it; Array::const_vector_iterator f_it; Array * filtered_J = NULL; if (filter_elements != empty_filter) { nb_element = filter_elements.getSize(); filtered_J = new Array(0, jac_loc.getNbComponent()); FEEngine::filterElementalData(mesh, jac_loc, *filtered_J, type, ghost_type, filter_elements); J_it = filtered_J->begin(); } else { nb_element = mesh.getNbElement(type, ghost_type); J_it = jac_loc.begin(); } intf.resize(nb_element * nb_points); f_it = in_f.begin(nb_degree_of_freedom); inte_it = intf.begin(nb_degree_of_freedom); for (UInt el = 0; el < nb_element; ++el, ++J_it, ++f_it, ++inte_it) { const Real & J = *J_it; const Vector & f = *f_it; Vector & inte_f = *inte_it; inte_f = f; inte_f *= J; } delete filtered_J; AKANTU_DEBUG_OUT(); } diff --git a/extra_packages/igfem/src/material_igfem/material_igfem_saw_tooth_damage.hh b/extra_packages/igfem/src/material_igfem/material_igfem_saw_tooth_damage.hh index 04ecc8063..479abdbd9 100644 --- a/extra_packages/igfem/src/material_igfem/material_igfem_saw_tooth_damage.hh +++ b/extra_packages/igfem/src/material_igfem/material_igfem_saw_tooth_damage.hh @@ -1,138 +1,138 @@ /** * @file material_igfem_saw_tooth_damage.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief Linear saw-tooth softening material model for IGFEM elements * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "material_damage.hh" #include "material_igfem_elastic.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MATERIAL_IGFEM_SAW_TOOTH_DAMAGE_HH_ #define AKANTU_MATERIAL_IGFEM_SAW_TOOTH_DAMAGE_HH_ namespace akantu { template class MaterialIGFEMSawToothDamage : public MaterialDamage { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ private: typedef MaterialDamage Parent; public: typedef std::pair ElementPair; MaterialIGFEMSawToothDamage(SolidMechanicsModel & model, const ID & id = ""); MaterialIGFEMSawToothDamage(SolidMechanicsModel & model, UInt dim, const Mesh & mesh, FEEngine & fe_engine, const ID & id = ""); virtual ~MaterialIGFEMSawToothDamage() {} protected: void initialize(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual void initMaterial(); /// virtual void updateInternalParameters(); virtual void computeAllStresses(GhostType ghost_type = _not_ghost); /// update internal field damage virtual UInt updateDamage(); - UInt updateDamage(UInt quad_index, const Real eq_stress, - ElementType el_type, GhostType ghost_type); + UInt updateDamage(UInt quad_index, const Real eq_stress, ElementType el_type, + GhostType ghost_type); /// update energies after damage has been updated // virtual void updateEnergiesAfterDamage(ElementType el_type, GhostType // ghost_typ); virtual void onBeginningSolveStep(const AnalysisMethod & method){}; virtual void onEndSolveStep(const AnalysisMethod & method){}; protected: /// constitutive law for all element of a type virtual void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); /// compute the equivalent stress on each Gauss point (i.e. the max prinicpal /// stress) and normalize it by the tensile strength virtual void computeNormalizedEquivalentStress(const Array & grad_u, ElementType el_type, GhostType ghost_type = _not_ghost); /// find max normalized equivalent stress void findMaxNormalizedEquivalentStress(ElementType el_type, GhostType ghost_type = _not_ghost); inline void computeDamageAndStressOnQuad(Matrix & sigma, Real & dam); protected: /* ------------------------------------------------------------------------ */ /* MeshEventHandler inherited members */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ virtual void onElementsAdded(const Array & element_list, const NewElementsEvent & event); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get max normalized equivalent stress AKANTU_GET_MACRO(NormMaxEquivalentStress, norm_max_equivalent_stress, Real); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// resistance to damage IGFEMInternalField Sc; /// internal field to store equivalent stress on each Gauss point IGFEMInternalField equivalent_stress; /// damage increment Real prescribed_dam; /// maximum equivalent stress Real norm_max_equivalent_stress; /// deviation from max stress at which Gauss point will still get damaged Real dam_tolerance; /// define damage threshold at which damage will be set to 1 Real dam_threshold; /// maximum damage value Real max_damage; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_igfem_saw_tooth_damage_inline_impl.hh" } // namespace akantu #endif /* AKANTU_MATERIAL_IGFEM_SAW_TOOTH_DAMAGE_HH_ */ diff --git a/extra_packages/igfem/src/shape_igfem.hh b/extra_packages/igfem/src/shape_igfem.hh index 537f6b401..93dd84432 100644 --- a/extra_packages/igfem/src/shape_igfem.hh +++ b/extra_packages/igfem/src/shape_igfem.hh @@ -1,199 +1,196 @@ /** * @file shape_igfem.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief shape functions for interface-enriched generalized FEM * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" #include "shape_functions.hh" #ifndef AKANTU_SHAPE_IGFEM_HH_ #define AKANTU_SHAPE_IGFEM_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ template <> class ShapeLagrange<_ek_igfem> : public ShapeFunctions { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: ShapeLagrange(const Mesh & mesh, const ID & id = "shape_igfem"); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: inline void initShapeFunctions(const Array & nodes, const Matrix & integration_points, const Matrix & integration_points_1, const Matrix & integration_points_2, - ElementType type, - GhostType ghost_type); + ElementType type, GhostType ghost_type); - inline void - interpolateEnrichmentsAllTypes(const Array & src, Array & dst, - ElementType type, - GhostType ghost_type) const; + inline void interpolateEnrichmentsAllTypes(const Array & src, + Array & dst, + ElementType type, + GhostType ghost_type) const; template inline void precomputeShapesOnEnrichedNodes(const Array & nodes, GhostType ghost_type); template void interpolateAtEnrichedNodes(const Array & src, Array & dst, GhostType ghost_type) const; /// pre compute all shapes on the element integration points from natural /// coordinates template void precomputeShapesOnIntegrationPoints(const Array & nodes, GhostType ghost_type); /// pre compute all shape derivatives on the element integration points from /// natural coordinates template void precomputeShapeDerivativesOnIntegrationPoints(const Array & nodes, GhostType ghost_type); /// interpolate nodal values on the integration points template void interpolateOnIntegrationPoints( const Array & u, Array & uq, UInt nb_degree_of_freedom, GhostType ghost_type = _not_ghost, const Array & filter_elements = empty_filter) const; /// interpolate on physical point template void interpolate(const Vector & real_coords, UInt elem, const Matrix & nodal_values, - Vector & interpolated, - GhostType ghost_type) const; + Vector & interpolated, GhostType ghost_type) const; /// compute the gradient of u on the integration points template void gradientOnIntegrationPoints( const Array & u, Array & nablauq, UInt nb_degree_of_freedom, GhostType ghost_type = _not_ghost, const Array & filter_elements = empty_filter) const; /// multiply a field by shape functions @f$ fts_{ij} = f_i * \varphi_j @f$ template void fieldTimesShapes(const Array & field, Array & field_times_shapes, GhostType ghost_type) const; /// find natural coords in parent element from real coords provided an element template void inverseMap(const Vector & real_coords, UInt element, Vector & natural_coords, GhostType ghost_type = _not_ghost) const; /// find natural coords in sub-element from real coords provided an element template void inverseMap(const Vector & real_coords, UInt element, Vector & natural_coords, UInt sub_element, GhostType ghost_type = _not_ghost) const; /// return true if the coordinates provided are inside the element, false /// otherwise template bool contains(const Vector & real_coords, UInt elem, GhostType ghost_type) const; /// compute the shape on a provided point template void computeShapes(const Vector & real_coords, UInt elem, Vector & shapes, GhostType ghost_type) const; /// compute the shape derivatives on a provided point template void computeShapeDerivatives(const Matrix & real_coords, UInt elem, Tensor3 & shapes, GhostType ghost_type) const; /// interpolate a field on a given physical point template void interpolateOnPhysicalPoint(const Vector & real_coords, UInt elem, const Array & field, Vector & interpolated, GhostType ghost_type) const; /// function to extract values at standard nodes and zero-out enriched values /// of a nodal field void extractValuesAtStandardNodes(const Array & nodal_values, Array & extracted_values, GhostType ghost_type) const; /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; protected: /// compute the shape derivatives on integration points for a given element template inline void computeShapeDerivativesOnCPointsByElement(const Matrix & node_coords, const Matrix & natural_coords, Tensor3 & shapesd) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get a the shapes vector - inline const Array & - getShapes(ElementType el_type, - GhostType ghost_type = _not_ghost) const; + inline const Array & getShapes(ElementType el_type, + GhostType ghost_type = _not_ghost) const; /// get a the shapes derivatives vector inline const Array & getShapesDerivatives(ElementType el_type, GhostType ghost_type = _not_ghost) const; /// get a the shapes vector inline const Array & getShapesAtEnrichedNodes(ElementType el_type, GhostType ghost_type = _not_ghost) const; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// shape functions for all elements ElementTypeMapArray shapes; /// shape functions derivatives for all elements ElementTypeMapArray shapes_derivatives; /// additional integration points for the IGFEM formulation ElementTypeMapArray igfem_integration_points; /// values of shape functions for all elements on the enriched nodes ElementTypeMapArray shapes_at_enrichments; }; } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "shape_igfem_inline_impl.hh" /// standard output stream operator // template // inline std::ostream & operator <<(std::ostream & stream, const // ShapeIGFEM & _this) // { // _this.printself(stream); // return stream; // } #endif /* AKANTU_SHAPE_IGFEM_HH_ */ diff --git a/extra_packages/igfem/src/shape_igfem_inline_impl.hh b/extra_packages/igfem/src/shape_igfem_inline_impl.hh index e72aaba33..a7febdfdf 100644 --- a/extra_packages/igfem/src/shape_igfem_inline_impl.hh +++ b/extra_packages/igfem/src/shape_igfem_inline_impl.hh @@ -1,733 +1,731 @@ /** * @file shape_igfem_inline_impl.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief ShapeIGFEM inline implementation * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SHAPE_IGFEM_INLINE_IMPL_HH_ #define AKANTU_SHAPE_IGFEM_INLINE_IMPL_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ inline const Array & ShapeLagrange<_ek_igfem>::getShapes(ElementType el_type, GhostType ghost_type) const { return shapes(FEEngine::getInterpolationType(el_type), ghost_type); } /* -------------------------------------------------------------------------- */ -inline const Array & ShapeLagrange<_ek_igfem>::getShapesDerivatives( - ElementType el_type, GhostType ghost_type) const { +inline const Array & +ShapeLagrange<_ek_igfem>::getShapesDerivatives(ElementType el_type, + GhostType ghost_type) const { return shapes_derivatives(FEEngine::getInterpolationType(el_type), ghost_type); } /* -------------------------------------------------------------------------- */ #define INIT_SHAPE_FUNCTIONS(type) \ setIntegrationPointsByType(integration_points, ghost_type); \ setIntegrationPointsByType::sub_element_type_1>( \ integration_points_1, ghost_type); \ setIntegrationPointsByType::sub_element_type_2>( \ integration_points_2, ghost_type); \ precomputeShapesOnIntegrationPoints(nodes, ghost_type); \ if (ElementClass::getNaturalSpaceDimension() == \ mesh.getSpatialDimension()) \ precomputeShapeDerivativesOnIntegrationPoints(nodes, ghost_type); \ precomputeShapesOnEnrichedNodes(nodes, ghost_type); inline void ShapeLagrange<_ek_igfem>::initShapeFunctions( const Array & nodes, const Matrix & integration_points, const Matrix & integration_points_1, const Matrix & integration_points_2, ElementType type, GhostType ghost_type) { AKANTU_BOOST_IGFEM_ELEMENT_SWITCH(INIT_SHAPE_FUNCTIONS); } #undef INIT_SHAPE_FUNCTIONS /* -------------------------------------------------------------------------- */ template inline void ShapeLagrange<_ek_igfem>::computeShapeDerivativesOnCPointsByElement( const Matrix & node_coords, const Matrix & natural_coords, Tensor3 & shapesd) const { AKANTU_DEBUG_IN(); // compute dnds Tensor3 dnds(node_coords.rows(), node_coords.cols(), natural_coords.cols()); ElementClass::computeDNDS(natural_coords, dnds); // compute dxds Tensor3 J(node_coords.rows(), natural_coords.rows(), natural_coords.cols()); ElementClass::computeJMat(dnds, node_coords, J); // compute shape derivatives ElementClass::computeShapeDerivatives(J, dnds, shapesd); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_igfem>::inverseMap(const Vector & real_coords, UInt elem, Vector & natural_coords, UInt sub_element, GhostType ghost_type) const { AKANTU_DEBUG_IN(); /// typedef for the two subelement_types and the parent element type const ElementType sub_type_1 = ElementClassProperty::sub_element_type_1; const ElementType sub_type_2 = ElementClassProperty::sub_element_type_2; UInt spatial_dimension = mesh.getSpatialDimension(); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); UInt * elem_val = mesh.getConnectivity(type, ghost_type).storage(); Matrix nodes_coord(spatial_dimension, nb_nodes_per_element); mesh.extractNodalValuesFromElement(mesh.getNodes(), nodes_coord.storage(), elem_val + elem * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); if (!sub_element) { UInt nb_nodes_sub_el = ElementClass::getNbNodesPerInterpolationElement(); Matrix sub_el_coords(spatial_dimension, nb_nodes_sub_el); ElementClass::getSubElementCoords(nodes_coord, sub_el_coords, sub_element); ElementClass::inverseMap(real_coords, sub_el_coords, natural_coords); } else { UInt nb_nodes_sub_el = ElementClass::getNbNodesPerInterpolationElement(); Matrix sub_el_coords(spatial_dimension, nb_nodes_sub_el); ElementClass::getSubElementCoords(nodes_coord, sub_el_coords, sub_element); ElementClass::inverseMap(real_coords, sub_el_coords, natural_coords); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_igfem>::inverseMap(const Vector & real_coords, UInt elem, Vector & natural_coords, GhostType ghost_type) const { /// map point into parent reference domain AKANTU_DEBUG_IN(); const ElementType parent_type = ElementClassProperty::parent_element_type; UInt spatial_dimension = mesh.getSpatialDimension(); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); UInt * elem_val = mesh.getConnectivity(type, ghost_type).storage(); Matrix nodes_coord(spatial_dimension, nb_nodes_per_element); mesh.extractNodalValuesFromElement(mesh.getNodes(), nodes_coord.storage(), elem_val + elem * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); UInt nb_nodes_parent_el = ElementClass::getNbNodesPerInterpolationElement(); Matrix parent_coords(spatial_dimension, nb_nodes_parent_el); ElementClass::getParentCoords(nodes_coord, parent_coords); ElementClass::inverseMap(real_coords, parent_coords, natural_coords); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template bool ShapeLagrange<_ek_igfem>::contains(const Vector & real_coords, - UInt elem, - GhostType ghost_type) const { + UInt elem, GhostType ghost_type) const { UInt spatial_dimension = mesh.getSpatialDimension(); Vector natural_coords(spatial_dimension); inverseMap(real_coords, elem, natural_coords, ghost_type); return ElementClass::contains(natural_coords); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_igfem>::interpolate(const Vector & real_coords, UInt elem, const Matrix & nodal_values, Vector & interpolated, GhostType ghost_type) const { UInt nb_shapes = ElementClass::getShapeSize(); Vector shapes(nb_shapes); computeShapes(real_coords, elem, shapes, ghost_type); ElementClass::interpolate(nodal_values, shapes, interpolated); } /* -------------------------------------------------------------------------- */ template -void ShapeLagrange<_ek_igfem>::computeShapes( - const Vector & real_coords, UInt elem, Vector & shapes, - GhostType ghost_type) const { +void ShapeLagrange<_ek_igfem>::computeShapes(const Vector & real_coords, + UInt elem, Vector & shapes, + GhostType ghost_type) const { AKANTU_DEBUG_IN(); /// typedef for the two subelement_types and the parent element type const ElementType sub_type_1 = ElementClassProperty::sub_element_type_1; const ElementType sub_type_2 = ElementClassProperty::sub_element_type_2; const ElementType parent_type = ElementClassProperty::parent_element_type; UInt spatial_dimension = mesh.getSpatialDimension(); /// parent contribution /// get the size of the parent shapes UInt size_of_parent_shapes = ElementClass::getShapeSize(); Vector parent_shapes(size_of_parent_shapes); /// compute parent shapes -> map shapes in the physical domain of the parent Vector natural_coords(spatial_dimension); Real tol = Math::getTolerance(); Math::setTolerance(1e-14); inverseMap(real_coords, elem, natural_coords, ghost_type); ElementClass::computeShapes(natural_coords, parent_shapes); natural_coords.zero(); /// sub-element contribution /// check which sub-element contains the physical point /// check if point is in sub-element 1 inverseMap(real_coords, elem, natural_coords, 0, ghost_type); if (ElementClass::contains(natural_coords)) { UInt size_of_sub_shapes = ElementClass::getShapeSize(); Vector sub_shapes(size_of_sub_shapes); ElementClass::computeShapes(natural_coords, sub_shapes); /// assemble shape functions ElementClass::assembleShapes(parent_shapes, sub_shapes, shapes, 0); } else { natural_coords.zero(); inverseMap(real_coords, elem, natural_coords, 1, ghost_type); AKANTU_DEBUG_ASSERT(ElementClass::contains(natural_coords), "Physical point not contained in any element"); UInt size_of_sub_shapes = ElementClass::getShapeSize(); Vector sub_shapes(size_of_sub_shapes); ElementClass::computeShapes(natural_coords, sub_shapes); /// assemble shape functions ElementClass::assembleShapes(parent_shapes, sub_shapes, shapes, 1); } Math::setTolerance(tol); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_igfem>::computeShapeDerivatives( const Matrix & real_coords, UInt elem, Tensor3 & shapesd, GhostType ghost_type) const { AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_igfem>::precomputeShapesOnIntegrationPoints( __attribute__((unused)) const Array & nodes, GhostType ghost_type) { AKANTU_DEBUG_IN(); InterpolationType itp_type = ElementClassProperty::interpolation_type; /// typedef for the two subelement_types and the parent element type const ElementType sub_type_1 = ElementClassProperty::sub_element_type_1; const ElementType sub_type_2 = ElementClassProperty::sub_element_type_2; const ElementType parent_type = ElementClassProperty::parent_element_type; /// get the spatial dimension for the given element type UInt spatial_dimension = ElementClass::getSpatialDimension(); /// get the integration points for the subelements Matrix & natural_coords_sub_1 = integration_points(sub_type_1, ghost_type); Matrix & natural_coords_sub_2 = integration_points(sub_type_2, ghost_type); /// store the number of quadrature points on each subelement and the toal /// number UInt nb_points_sub_1 = natural_coords_sub_1.cols(); UInt nb_points_sub_2 = natural_coords_sub_2.cols(); UInt nb_total_points = nb_points_sub_1 + nb_points_sub_2; // get the integration points for the parent element UInt nb_element = mesh.getConnectivity(type, ghost_type).getSize(); Array & natural_coords_parent = igfem_integration_points.alloc( nb_element * nb_total_points, spatial_dimension, type, ghost_type); Array::matrix_iterator natural_coords_parent_it = natural_coords_parent.begin_reinterpret(spatial_dimension, nb_total_points, nb_element); /// get the size of the shapes UInt size_of_shapes = ElementClass::getShapeSize(); UInt size_of_parent_shapes = ElementClass::getShapeSize(); UInt size_of_sub_1_shapes = ElementClass::getShapeSize(); UInt size_of_sub_2_shapes = ElementClass::getShapeSize(); /// initialize the matrices to store the shape functions of the subelements /// and the parent Matrix sub_1_shapes(size_of_sub_1_shapes, nb_points_sub_1); Matrix sub_2_shapes(size_of_sub_2_shapes, nb_points_sub_2); Matrix parent_1_shapes(size_of_parent_shapes, nb_points_sub_1); Matrix parent_2_shapes(size_of_parent_shapes, nb_points_sub_2); /// compute the shape functions of the subelements ElementClass::computeShapes(natural_coords_sub_1, sub_1_shapes); ElementClass::computeShapes(natural_coords_sub_2, sub_2_shapes); /// get the nodal coordinates per element UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); Array x_el(0, spatial_dimension * nb_nodes_per_element); FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type); Array::matrix_iterator x_it = x_el.begin(spatial_dimension, nb_nodes_per_element); /// allocate the shapes for the given element type Array & shapes_tmp = shapes.alloc(nb_element * nb_total_points, size_of_shapes, itp_type, ghost_type); Array::matrix_iterator shapes_it = shapes_tmp.begin_reinterpret( ElementClass::getNbNodesPerInterpolationElement(), nb_total_points, nb_element); Matrix physical_points_1(spatial_dimension, nb_points_sub_1); Matrix physical_points_2(spatial_dimension, nb_points_sub_2); Matrix parent_natural_coords_1(spatial_dimension, nb_points_sub_1); Matrix parent_natural_coords_2(spatial_dimension, nb_points_sub_2); /// intialize the matrices for the parent and subelement coordinates UInt nb_nodes_parent_el = ElementClass::getNbNodesPerInterpolationElement(); UInt nb_nodes_sub_el_1 = ElementClass::getNbNodesPerInterpolationElement(); UInt nb_nodes_sub_el_2 = ElementClass::getNbNodesPerInterpolationElement(); Matrix parent_coords(spatial_dimension, nb_nodes_parent_el); Matrix sub_el_1_coords(spatial_dimension, nb_nodes_sub_el_1); Matrix sub_el_2_coords(spatial_dimension, nb_nodes_sub_el_2); /// loop over all elements of the given type and compute the shape functions Vector all_shapes(size_of_shapes); for (UInt elem = 0; elem < nb_element; ++elem, ++shapes_it, ++x_it, ++natural_coords_parent_it) { Matrix & N = *shapes_it; const Matrix & X = *x_it; Matrix & nc_parent = *natural_coords_parent_it; /// map the sub element integration points into the parent reference domain ElementClass::mapFromSubRefToParentRef( X, sub_el_1_coords, parent_coords, sub_1_shapes, physical_points_1, parent_natural_coords_1, 0); ElementClass::mapFromSubRefToParentRef( X, sub_el_2_coords, parent_coords, sub_2_shapes, physical_points_2, parent_natural_coords_2, 1); /// compute the parent shape functions on all integration points ElementClass::computeShapes(parent_natural_coords_1, parent_1_shapes); ElementClass::computeShapes(parent_natural_coords_2, parent_2_shapes); /// copy the results into the shape functions iterator and natural coords /// iterator for (UInt i = 0; i < nb_points_sub_1; ++i) { ElementClass::assembleShapes(parent_1_shapes(i), sub_1_shapes(i), all_shapes, 0); N(i) = all_shapes; nc_parent(i) = parent_natural_coords_1(i); } for (UInt i = 0; i < nb_points_sub_2; ++i) { ElementClass::assembleShapes(parent_2_shapes(i), sub_2_shapes(i), all_shapes, 1); N(i + nb_points_sub_1) = all_shapes; /// N(i + nb_points_sub_2) = all_shapes; nc_parent(i + nb_points_sub_1) = parent_natural_coords_2(i); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_igfem>::precomputeShapeDerivativesOnIntegrationPoints( const Array & nodes, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// typedef for the two subelement_types and the parent element type const ElementType sub_type_1 = ElementClassProperty::sub_element_type_1; const ElementType sub_type_2 = ElementClassProperty::sub_element_type_2; const ElementType parent_type = ElementClassProperty::parent_element_type; InterpolationType itp_type = ElementClassProperty::interpolation_type; UInt spatial_dimension = mesh.getSpatialDimension(); /// get the integration points for the subelements Matrix & natural_coords_sub_1 = integration_points(sub_type_1, ghost_type); Matrix & natural_coords_sub_2 = integration_points(sub_type_2, ghost_type); /// store the number of quadrature points on each subelement and the toal /// number UInt nb_points_sub_1 = natural_coords_sub_1.cols(); UInt nb_points_sub_2 = natural_coords_sub_2.cols(); UInt nb_points_total = nb_points_sub_1 + nb_points_sub_2; UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); UInt size_of_shapesd = ElementClass::getShapeDerivativesSize(); /// intialize the matrices for the parent and subelement coordinates UInt nb_nodes_parent_el = ElementClass::getNbNodesPerInterpolationElement(); UInt nb_nodes_sub_el_1 = ElementClass::getNbNodesPerInterpolationElement(); UInt nb_nodes_sub_el_2 = ElementClass::getNbNodesPerInterpolationElement(); Matrix parent_coords(spatial_dimension, nb_nodes_parent_el); Matrix sub_el_1_coords(spatial_dimension, nb_nodes_sub_el_1); Matrix sub_el_2_coords(spatial_dimension, nb_nodes_sub_el_2); UInt nb_element = mesh.getConnectivity(type, ghost_type).getSize(); Array & shapes_derivatives_tmp = shapes_derivatives.alloc( nb_element * nb_points_total, size_of_shapesd, itp_type, ghost_type); /// get an iterator to the coordiantes of the elements Array x_el(0, spatial_dimension * nb_nodes_per_element); FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type); Real * shapesd_val = shapes_derivatives_tmp.storage(); Array::matrix_iterator x_it = x_el.begin(spatial_dimension, nb_nodes_per_element); /// get an iterator to the integration points of the parent element Array & natural_coords_parent = igfem_integration_points(type, ghost_type); Array::matrix_iterator natural_coords_parent_it = natural_coords_parent.begin_reinterpret(spatial_dimension, nb_points_total, nb_element); Tensor3 B_sub_1(spatial_dimension, nb_nodes_sub_el_1, nb_points_sub_1); Tensor3 B_sub_2(spatial_dimension, nb_nodes_sub_el_2, nb_points_sub_2); Tensor3 B_parent(spatial_dimension, nb_nodes_parent_el, nb_points_total); /// assemble the shape derivatives Matrix all_shapes(spatial_dimension, nb_nodes_per_element); for (UInt elem = 0; elem < nb_element; ++elem, ++x_it, ++natural_coords_parent_it) { Matrix & X = *x_it; Matrix & nc_parent = *natural_coords_parent_it; Tensor3 B(shapesd_val, spatial_dimension, nb_nodes_per_element, nb_points_total); /// get the coordinates of the two sub elements and the parent element ElementClass::getSubElementCoords(X, sub_el_1_coords, 0); ElementClass::getSubElementCoords(X, sub_el_2_coords, 1); ElementClass::getParentCoords(X, parent_coords); /// compute the subelements' shape derivatives and the parent shape /// derivatives computeShapeDerivativesOnCPointsByElement( sub_el_1_coords, natural_coords_sub_1, B_sub_1); computeShapeDerivativesOnCPointsByElement( sub_el_2_coords, natural_coords_sub_2, B_sub_2); computeShapeDerivativesOnCPointsByElement(parent_coords, nc_parent, B_parent); for (UInt i = 0; i < nb_points_sub_1; ++i) { ElementClass::assembleShapeDerivatives(B_parent(i), B_sub_1(i), all_shapes, 0); B(i) = all_shapes; } for (UInt i = 0; i < nb_points_sub_2; ++i) { ElementClass::assembleShapeDerivatives(B_parent(i), B_sub_2(i), all_shapes, 1); B(i + nb_points_sub_1) = all_shapes; } shapesd_val += size_of_shapesd * nb_points_total; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_igfem>::interpolateOnIntegrationPoints( const Array & in_u, Array & out_uq, UInt nb_degree_of_freedom, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); InterpolationType itp_type = ElementClassProperty::interpolation_type; AKANTU_DEBUG_ASSERT(shapes.exists(itp_type, ghost_type), "No shapes for the type " << shapes.printType(itp_type, ghost_type)); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); Array u_el(0, nb_degree_of_freedom * nb_nodes_per_element); FEEngine::extractNodalToElementField(mesh, in_u, u_el, type, ghost_type, filter_elements); this->interpolateElementalFieldOnIntegrationPoints( u_el, out_uq, ghost_type, shapes(itp_type, ghost_type), filter_elements); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_igfem>::gradientOnIntegrationPoints( const Array & in_u, Array & out_nablauq, UInt nb_degree_of_freedom, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); InterpolationType itp_type = ElementClassProperty::interpolation_type; AKANTU_DEBUG_ASSERT( shapes_derivatives.exists(itp_type, ghost_type), "No shapes derivatives for the type " << shapes_derivatives.printType(itp_type, ghost_type)); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); Array u_el(0, nb_degree_of_freedom * nb_nodes_per_element); FEEngine::extractNodalToElementField(mesh, in_u, u_el, type, ghost_type, filter_elements); this->gradientElementalFieldOnIntegrationPoints( u_el, out_nablauq, ghost_type, shapes_derivatives(itp_type, ghost_type), filter_elements); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_igfem>::fieldTimesShapes( const Array & field, Array & field_times_shapes, GhostType ghost_type) const { AKANTU_DEBUG_IN(); field_times_shapes.resize(field.getSize()); UInt size_of_shapes = ElementClass::getShapeSize(); InterpolationType itp_type = ElementClassProperty::interpolation_type; UInt nb_degree_of_freedom = field.getNbComponent(); const Array & shape = shapes(itp_type, ghost_type); Array::const_matrix_iterator field_it = field.begin(nb_degree_of_freedom, 1); Array::const_matrix_iterator shapes_it = shape.begin(1, size_of_shapes); Array::matrix_iterator it = field_times_shapes.begin(nb_degree_of_freedom, size_of_shapes); Array::matrix_iterator end = field_times_shapes.end(nb_degree_of_freedom, size_of_shapes); for (; it != end; ++it, ++field_it, ++shapes_it) { it->mul(*field_it, *shapes_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_igfem>::interpolateOnPhysicalPoint( const Vector & real_coords, UInt elem, const Array & field, Vector & interpolated, GhostType ghost_type) const { AKANTU_DEBUG_IN(); Vector shapes(ElementClass::getShapeSize()); computeShapes(real_coords, elem, shapes, ghost_type); UInt spatial_dimension = mesh.getSpatialDimension(); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); UInt * elem_val = mesh.getConnectivity(type, ghost_type).storage(); Matrix nodes_val(spatial_dimension, nb_nodes_per_element); mesh.extractNodalValuesFromElement(field, nodes_val.storage(), elem_val + elem * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); ElementClass::interpolate(nodes_val, shapes, interpolated); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_igfem>::precomputeShapesOnEnrichedNodes( - __attribute__((unused)) const Array & nodes, - GhostType ghost_type) { + __attribute__((unused)) const Array & nodes, GhostType ghost_type) { AKANTU_DEBUG_IN(); InterpolationType itp_type = ElementClassProperty::interpolation_type; const ElementType parent_type = ElementClassProperty::parent_element_type; const ElementType sub_type = ElementClassProperty::sub_element_type_1; /// get the spatial dimension for the given element type UInt spatial_dimension = ElementClass::getSpatialDimension(); // get the integration points for the parent element UInt nb_element = mesh.getConnectivity(type, ghost_type).getSize(); /// get the size of the shapes UInt nb_enriched_nodes = ElementClass::getNbEnrichments(); UInt nb_parent_nodes = ElementClass::getNbNodesPerInterpolationElement(); UInt size_of_shapes = ElementClass::getShapeSize(); UInt size_of_parent_shapes = ElementClass::getShapeSize(); UInt size_of_sub_shapes = ElementClass::getShapeSize(); Vector parent_shapes(size_of_parent_shapes); Vector sub_shapes(size_of_sub_shapes); Vector shapes(size_of_shapes); /// get the nodal coordinates per element UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); Array x_el(0, spatial_dimension * nb_nodes_per_element); FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type); Array::matrix_iterator x_it = x_el.begin(spatial_dimension, nb_nodes_per_element); /// allocate the shapes for the given element type Array & shapes_tmp = shapes_at_enrichments.alloc( nb_element * nb_enriched_nodes, size_of_shapes, itp_type, ghost_type); Array::matrix_iterator shapes_it = shapes_tmp.begin_reinterpret( ElementClass::getNbNodesPerInterpolationElement(), nb_enriched_nodes, nb_element); Vector real_coords(spatial_dimension); Vector natural_coords(spatial_dimension); Matrix parent_coords(spatial_dimension, nb_parent_nodes); UInt * sub_element_enrichments = ElementClass::getSubElementEnrichments(); /// loop over all elements for (UInt elem = 0; elem < nb_element; ++elem, ++shapes_it, ++x_it) { Matrix & N = *shapes_it; const Matrix & X = *x_it; for (UInt i = 0; i < nb_enriched_nodes; ++i) { /// get the parent element coordinates ElementClass::getParentCoords(X, parent_coords); /// get the physical coords of the enriched node real_coords = X(nb_parent_nodes + i); /// map the physical point into the parent ref domain ElementClass::inverseMap(real_coords, parent_coords, natural_coords); /// compute the parent shape functions ElementClass::computeShapes(natural_coords, parent_shapes); /// Sub-element contribution sub_shapes.zero(); sub_shapes(sub_element_enrichments[i]) = 1.; ElementClass::assembleShapes(parent_shapes, sub_shapes, shapes, 0); N(i) = shapes; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_igfem>::interpolateAtEnrichedNodes( - const Array & src, Array & dst, - GhostType ghost_type) const { + const Array & src, Array & dst, GhostType ghost_type) const { AKANTU_DEBUG_IN(); const ElementType parent_type = ElementClassProperty::parent_element_type; UInt nb_element = mesh.getNbElement(type, ghost_type); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); UInt nb_parent_nodes = ElementClass::getNbNodesPerInterpolationElement(); UInt nb_enrichments = ElementClass::getNbEnrichments(); UInt * elem_val = mesh.getConnectivity(type, ghost_type).storage(); UInt spatial_dimension = mesh.getSpatialDimension(); Matrix nodes_val(spatial_dimension, nb_nodes_per_element); InterpolationType itp_type = ElementClassProperty::interpolation_type; const Array & shapes = shapes_at_enrichments(itp_type, ghost_type); Array::const_matrix_iterator shapes_it = shapes.begin_reinterpret( nb_nodes_per_element, nb_enrichments, nb_element); Array::vector_iterator dst_vect = dst.begin(spatial_dimension); Vector interpolated(spatial_dimension); for (UInt e = 0; e < nb_element; ++e, ++shapes_it) { const Matrix & el_shapes = *shapes_it; mesh.extractNodalValuesFromElement(src, nodes_val.storage(), elem_val + e * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); ; for (UInt i = 0; i < nb_enrichments; ++i) { ElementClass::interpolate(nodes_val, el_shapes(i), interpolated); UInt enr_node_idx = elem_val[e * nb_nodes_per_element + nb_parent_nodes + i]; dst_vect[enr_node_idx] = interpolated; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ #define COMPUTE_ENRICHED_VALUES(type) \ interpolateAtEnrichedNodes(src, dst, ghost_type); inline void ShapeLagrange<_ek_igfem>::interpolateEnrichmentsAllTypes( const Array & src, Array & dst, ElementType type, GhostType ghost_type) const { AKANTU_BOOST_IGFEM_ELEMENT_SWITCH(COMPUTE_ENRICHED_VALUES); } #undef COMPUTE_ENRICHED_VALUES /* -------------------------------------------------------------------------- */ } // namespace akantu #endif /* AKANTU_SHAPE_IGFEM_INLINE_IMPL_HH_ */ diff --git a/extra_packages/igfem/src/solid_mechanics_model_igfem.hh b/extra_packages/igfem/src/solid_mechanics_model_igfem.hh index 4933932f8..c9f893958 100644 --- a/extra_packages/igfem/src/solid_mechanics_model_igfem.hh +++ b/extra_packages/igfem/src/solid_mechanics_model_igfem.hh @@ -1,196 +1,196 @@ /** * @file solid_mechanics_model_igfem.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief solid mechanics model for IGFEM analysis * * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SOLID_MECHANICS_MODEL_IGFEM_HH_ #define AKANTU_SOLID_MECHANICS_MODEL_IGFEM_HH_ #include "global_ids_updater.hh" #include "igfem_enrichment.hh" #include "solid_mechanics_model.hh" #include "solid_mechanics_model_event_handler.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ struct SolidMechanicsModelIGFEMOptions : public SolidMechanicsModelOptions { SolidMechanicsModelIGFEMOptions(AnalysisMethod analysis_method = _static, bool no_init_materials = false) : SolidMechanicsModelOptions(analysis_method, no_init_materials) {} }; extern const SolidMechanicsModelIGFEMOptions default_solid_mechanics_model_igfem_options; /* -------------------------------------------------------------------------- */ /* Solid Mechanics Model for IGFEM analysis */ /* -------------------------------------------------------------------------- */ class SolidMechanicsModelIGFEM : public SolidMechanicsModel, public SolidMechanicsModelEventHandler, public IGFEMEnrichment { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef FEEngineTemplate MyFEEngineIGFEMType; typedef std::map ElementMap; typedef std::map FEEnginesPerKindMap; SolidMechanicsModelIGFEM(Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "solid_mechanics_model_igfem"); virtual ~SolidMechanicsModelIGFEM(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize the cohesive model virtual void initFull(const ModelOptions & options = default_solid_mechanics_model_igfem_options); /// initialize the model virtual void initModel(); /// initialize igfem material virtual void initMaterials(); /// register the tags associated with the parallel synchronizer virtual void initParallel(MeshPartition * partition, DataAccessor * data_accessor = NULL); /// allocate all vectors virtual void initArrays(); /// transfer internals from old to new elements void transferInternalValues(const ID & internal, std::vector & new_elements, Array & added_quads, Array & internal_values); /// compute the barycenter for a sub-element inline void getSubElementBarycenter(UInt element, UInt sub_element, ElementType type, Vector & barycenter, GhostType ghost_type) const; /// apply a constant eigen_grad_u on all quadrature points of a given material virtual void applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const ID & material_name, const GhostType ghost_type = _not_ghost); private: /// compute the real values of displacement, force, etc. on the enriched nodes void computeValuesOnEnrichedNodes(); /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: virtual void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event); virtual void onNodesRemoved(const Array & element_list, const Array & new_numbering, const RemovedNodesEvent & event); virtual void onElementsAdded(const Array & nodes_list, const NewElementsEvent & event); virtual void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event); /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: virtual void onDump(); virtual void addDumpGroupFieldToDumper(const std::string & dumper_name, const std::string & field_id, const std::string & group_name, ElementKind element_kind, bool padding_flag); virtual dumpers::Field * createElementalField(const std::string & field_name, - const std::string & group_name, - bool padding_flag, - const UInt & spatial_dimension, - ElementKind kind); + const std::string & group_name, + bool padding_flag, + const UInt & spatial_dimension, + ElementKind kind); virtual dumpers::Field * createNodalFieldReal(const std::string & field_name, - const std::string & group_name, - bool padding_flag); + const std::string & group_name, + bool padding_flag); /* -------------------------------------------------------------------------- */ /* Accessors */ /* -------------------------------------------------------------------------- */ public: /// get the fe-engines per kind AKANTU_GET_MACRO(FEEnginesPerKind, fe_engines_per_kind, const FEEnginesPerKindMap &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// real displacements array Array * real_displacement; /// real forces array Array * real_force; /// real residuals array Array * real_residual; /// map between and new elements (needed when the interface is moving) ElementMap element_map; /// global connectivity ids updater GlobalIdsUpdater * global_ids_updater; /// map between element kind and corresponding FEEngine object FEEnginesPerKindMap fe_engines_per_kind; }; /* -------------------------------------------------------------------------- */ /* IGFEMMaterialSelector */ /* -------------------------------------------------------------------------- */ class DefaultMaterialIGFEMSelector : public DefaultMaterialSelector { public: DefaultMaterialIGFEMSelector(const SolidMechanicsModelIGFEM & model) : DefaultMaterialSelector(model.getMaterialByElement()), fallback_value_igfem(0) {} virtual UInt operator()(const Element & element) { if (Mesh::getKind(element.type) == _ek_igfem) return fallback_value_igfem; else return DefaultMaterialSelector::operator()(element); } void setIGFEMFallback(UInt f) { this->fallback_value_igfem = f; } protected: UInt fallback_value_igfem; }; } // namespace akantu #if defined(AKANTU_INCLUDE_INLINE_IMPL) #include "solid_mechanics_model_igfem_inline_impl.hh" #endif #endif /* AKANTU_SOLID_MECHANICS_MODEL_IGFEM_HH_ */ diff --git a/extra_packages/igfem/src/solid_mechanics_model_igfem_inline_impl.hh b/extra_packages/igfem/src/solid_mechanics_model_igfem_inline_impl.hh index 4c708c938..3158409f9 100644 --- a/extra_packages/igfem/src/solid_mechanics_model_igfem_inline_impl.hh +++ b/extra_packages/igfem/src/solid_mechanics_model_igfem_inline_impl.hh @@ -1,59 +1,59 @@ /** * @file solid_mechanics_model_igfem_inline_impl.hh * @author Aurelia Isabel Cuba Ramos * @date Wed Nov 4 15:53:52 2015 * * @brief Implementation on inline functions for SMMIGFEM * * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 of the License, or (at your option) any * later version. * * Akantu 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 Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SOLID_MECHANICS_MODEL_IGFEM_INLINE_IMPL_HH_ #define AKANTU_SOLID_MECHANICS_MODEL_IGFEM_INLINE_IMPL_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ inline void SolidMechanicsModelIGFEM::getSubElementBarycenter( - UInt element, UInt sub_element, ElementType type, - Vector & barycenter, GhostType ghost_type) const { + UInt element, UInt sub_element, ElementType type, Vector & barycenter, + GhostType ghost_type) const { UInt * conn_val = this->mesh.getConnectivity(type, ghost_type).storage(); UInt nb_sub_element_nodes = IGFEMHelper::getNbNodesPerSubElement(type, sub_element); UInt * sub_el_conn = IGFEMHelper::getSubElementConnectivity(type, sub_element); UInt nb_nodes_per_element = this->mesh.getNbNodesPerElement(type); const Array & node_coords = this->mesh.getNodes(); Real local_coord[spatial_dimension * nb_sub_element_nodes]; UInt offset = element * nb_nodes_per_element; for (UInt n = 0; n < nb_sub_element_nodes; ++n) { UInt index = conn_val[offset + sub_el_conn[n]]; memcpy(local_coord + n * spatial_dimension, node_coords.storage() + index * spatial_dimension, spatial_dimension * sizeof(Real)); } Math::barycenter(local_coord, nb_sub_element_nodes, spatial_dimension, barycenter.storage()); } } // namespace akantu #endif /* AKANTU_SOLID_MECHANICS_MODEL_IGFEM_INLINE_IMPL_HH_ */ diff --git a/extra_packages/traction-at-split-node-contact/src/boundary_conditions/force_based_dirichlet.hh b/extra_packages/traction-at-split-node-contact/src/boundary_conditions/force_based_dirichlet.hh index 3839869bf..a859105fa 100644 --- a/extra_packages/traction-at-split-node-contact/src/boundary_conditions/force_based_dirichlet.hh +++ b/extra_packages/traction-at-split-node-contact/src/boundary_conditions/force_based_dirichlet.hh @@ -1,130 +1,130 @@ /** * @file force_based_dirichlet.hh * * @author Dana Christen * @author David Simon Kammer * * @date creation: Fri Mar 16 2018 * @date last modification: Tue Sep 29 2020 * * @brief dirichlet boundary condition that tries to keep the force at a given * value * * * @section LICENSE * * Copyright (©) 2015-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 of the License, or (at your option) any * later version. - * + * * Akantu 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 Lesser General Public License for more * details. - * + * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef AST_FORCE_BASED_DIRICHLET_HH_ #define AST_FORCE_BASED_DIRICHLET_HH_ // akantu #include "aka_common.hh" namespace akantu { /* -------------------------------------------------------------------------- */ class ForceBasedDirichlet : public BC::Dirichlet::IncrementValue { protected: typedef const Array * RealArrayPtr; typedef const Array * IntArrayPtr; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: ForceBasedDirichlet(SolidMechanicsModel & model, BC::Axis ax, Real target_f, Real mass = 0.) : IncrementValue(0., ax), model(model), mass(mass), velocity(0.), target_force(target_f), total_residual(0.) {} virtual ~ForceBasedDirichlet() {} /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void updateTotalResidual() { this->total_residual = 0.; for (auto && subboundary : this->subboundaries) { this->total_residual += integrateResidual(subboundary, this->model, this->axis); } } virtual Real update() { AKANTU_DEBUG_IN(); this->updateTotalResidual(); Real total_force = this->target_force + this->total_residual; Real a = total_force / this->mass; Real dt = model.getTimeStep(); this->velocity += 0.5 * dt * a; this->value = this->velocity * dt + 0.5 * dt * dt * a; // increment position dx this->velocity += 0.5 * dt * a; AKANTU_DEBUG_OUT(); return this->total_residual; } Real applyYourself() { AKANTU_DEBUG_IN(); Real reaction = this->update(); for (auto && subboundary : this->subboundaries) { this->model.applyBC(*this, subboundary); } AKANTU_DEBUG_OUT(); return reaction; } /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_SET_MACRO(Mass, mass, Real); AKANTU_SET_MACRO(TargetForce, target_force, Real); void insertSubBoundary(const std::string & sb_name) { this->subboundaries.insert(sb_name); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ typedef std::set SubBoundarySet; protected: SolidMechanicsModel & model; SubBoundarySet subboundaries; Real mass; Real velocity; Real target_force; Real total_residual; }; } // namespace akantu #endif /* AST_FORCE_BASED_DIRICHLET_HH_ */ diff --git a/extra_packages/traction-at-split-node-contact/src/boundary_conditions/inclined_flat_dirichlet.hh b/extra_packages/traction-at-split-node-contact/src/boundary_conditions/inclined_flat_dirichlet.hh index 513ea3549..10b378ee2 100644 --- a/extra_packages/traction-at-split-node-contact/src/boundary_conditions/inclined_flat_dirichlet.hh +++ b/extra_packages/traction-at-split-node-contact/src/boundary_conditions/inclined_flat_dirichlet.hh @@ -1,82 +1,82 @@ /** * @file inclined_flat_dirichlet.hh * * @author David Simon Kammer * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Sep 29 2020 * * @brief inclined dirichlet * * * @section LICENSE * * Copyright (©) 2015-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 of the License, or (at your option) any * later version. - * + * * Akantu 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 Lesser General Public License for more * details. - * + * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef AST_INCLINED_FLAT_DIRICHLET_HH_ #define AST_INCLINED_FLAT_DIRICHLET_HH_ // akantu #include "aka_common.hh" namespace akantu { /* -------------------------------------------------------------------------- */ class InclinedFlatDirichlet : public BC::Dirichlet::DirichletFunctor { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: InclinedFlatDirichlet(Real val, BC::Axis ax, BC::Axis incl_ax, Real center_coord, Real tang) : DirichletFunctor(ax), value(val), incl_ax(incl_ax), center_coord(center_coord), tang(tang){}; virtual ~InclinedFlatDirichlet() {} /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: inline void operator()(UInt node, Vector & flags, Vector & primal, const Vector & coord) const { AKANTU_DEBUG_IN(); Real dist = coord(incl_ax) - this->center_coord; flags(axis) = true; primal(axis) = this->value + this->tang * dist; AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: Real value; BC::Axis incl_ax; Real center_coord; Real tang; }; } // namespace akantu #endif /* AST_INCLINED_FLAT_DIRICHLET_HH_ */ diff --git a/extra_packages/traction-at-split-node-contact/src/boundary_conditions/spring_bc.hh b/extra_packages/traction-at-split-node-contact/src/boundary_conditions/spring_bc.hh index e6105c662..05afbadbc 100644 --- a/extra_packages/traction-at-split-node-contact/src/boundary_conditions/spring_bc.hh +++ b/extra_packages/traction-at-split-node-contact/src/boundary_conditions/spring_bc.hh @@ -1,143 +1,143 @@ /** * @file spring_bc.hh * * @author Dana Christen * @author David Simon Kammer * * @date creation: Fri Mar 16 2018 * @date last modification: Tue Sep 29 2020 * * @brief spring boundary condition * * * @section LICENSE * * Copyright (©) 2015-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 of the License, or (at your option) any * later version. - * + * * Akantu 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 Lesser General Public License for more * details. - * + * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef AST_SPRING_BC_HH_ #define AST_SPRING_BC_HH_ // simtools #include "force_based_dirichlet.hh" namespace akantu { /* -------------------------------------------------------------------------- */ class SpringBC : public ForceBasedDirichlet { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: SpringBC(SolidMechanicsModel & model, BC::Axis ax, Real stiffness, Real mass = 0.) : ForceBasedDirichlet(model, ax, 0., mass), stiffness(stiffness), elongation(0.) {} virtual ~SpringBC() {} /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual Real update() { AKANTU_DEBUG_IN(); this->target_force = -this->stiffness * this->elongation; Real reaction = ForceBasedDirichlet::update(); this->elongation += this->value; AKANTU_DEBUG_OUT(); return reaction; } /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Elongation, elongation, Real); inline void setToEquilibrium() { AKANTU_DEBUG_IN(); this->updateTotalResidual(); this->target_force = -this->total_residual; this->elongation = -this->target_force / this->stiffness; AKANTU_DEBUG_OUT(); } /// change elongation /// dx > 0 -> target_force < 0 inline void incrementElongation(Real dx) { AKANTU_DEBUG_IN(); this->elongation += dx; AKANTU_DEBUG_OUT(); } // friend std::ostream& operator<<(std::ostream& out, const SpringBC & // spring); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: Real stiffness; Real elongation; }; // class SpringBCRestricted : public SpringBC { // public: // SpringBCRestricted(BC::Axis ax, Real target_force, BC::Axis surface_axis, // Real min, Real max) // :SpringBC(ax, target_force), surface_axis(surface_axis), min(min), // max(max) {} // virtual ~SpringBCRestricted() {} // public: // inline void operator()(UInt node, Vector & flags, Vector & // primal, const Vector & coord) const { // if(coord(surface_axis) > min && coord(surface_axis) < max) { // SpringBC::operator()(node, flags, primal, coord); // } // } // private: // BC::Axis surface_axis; // Real min; // Real max; // }; // std::ostream& operator<<(std::ostream& out, const SpringBC & spring) { // out << "Real total_residual: " << *spring.total_residual << std::endl; // out << "Real mass: " << spring.mass << std::endl; // out << "Real k: " << spring.k << std::endl; // out << "Real delta_x: " << spring.delta_x << std::endl; // out << "Real dt: " << spring.dt << std::endl; // out << "Real v: " << spring.v << std::endl; // out << "Real dx: " << spring.dx << std::endl; // out << "Real forcing_vel: " << spring.forcing_vel << std::endl; // return out; // } } // namespace akantu #endif /* AST_SPRING_BC_HH_ */ diff --git a/extra_packages/traction-at-split-node-contact/src/common/manual_restart.hh b/extra_packages/traction-at-split-node-contact/src/common/manual_restart.hh index 5ee0ff94e..26d7126e0 100644 --- a/extra_packages/traction-at-split-node-contact/src/common/manual_restart.hh +++ b/extra_packages/traction-at-split-node-contact/src/common/manual_restart.hh @@ -1,53 +1,53 @@ /** * @file manual_restart.hh * * @author David Simon Kammer * * @date creation: Fri Mar 16 2018 * @date last modification: Fri Mar 16 2018 * * @brief Restart tools header * * * @section LICENSE * * Copyright (©) 2015-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 of the License, or (at your option) any * later version. - * + * * Akantu 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 Lesser General Public License for more * details. - * + * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /** * @file manual_restart.hh * @author Dana Christen * @date May 15, 2013 */ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" #include "solid_mechanics_model.hh" void dumpArray(const akantu::Array & array, const std::string & fname); void loadArray(akantu::Array & array, const std::string & fname); void loadRestart(akantu::SolidMechanicsModel & model, const std::string & fname, akantu::UInt prank); void loadRestart(akantu::SolidMechanicsModel & model, const std::string & fname); void dumpRestart(akantu::SolidMechanicsModel & model, const std::string & fname, akantu::UInt prank); void dumpRestart(akantu::SolidMechanicsModel & model, const std::string & fname); diff --git a/extra_packages/traction-at-split-node-contact/src/common/parameter_reader.hh b/extra_packages/traction-at-split-node-contact/src/common/parameter_reader.hh index 321fe2523..d8e9ef118 100644 --- a/extra_packages/traction-at-split-node-contact/src/common/parameter_reader.hh +++ b/extra_packages/traction-at-split-node-contact/src/common/parameter_reader.hh @@ -1,111 +1,111 @@ /** * @file parameter_reader.hh * * @author David Simon Kammer * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Sep 29 2020 * * @brief for simulations to read parameters from an input file * * * @section LICENSE * * Copyright (©) 2015-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 of the License, or (at your option) any * later version. - * + * * Akantu 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 Lesser General Public License for more * details. - * + * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef AST_PARAMETER_READER_HH_ #define AST_PARAMETER_READER_HH_ /* -------------------------------------------------------------------------- */ // std #include #include // akantu #include "aka_common.hh" namespace akantu { /* -------------------------------------------------------------------------- */ class ParameterReader { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: ParameterReader(); virtual ~ParameterReader() = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// read input file void readInputFile(std::string file_name); /// write input file void writeInputFile(std::string file_name) const; /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// template T get(std::string key) const; template bool has(std::string key) const; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// type of data available std::set data_types; /// data std::map element_type_data; std::map string_data; std::map int_data; std::map uint_data; std::map real_data; std::map bool_data; /// convert string to element type std::map _input_to_akantu_element_types; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "parameter_reader_inline_impl.hh" /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const ParameterReader & _this) { _this.printself(stream); return stream; } } // namespace akantu #endif /* AST_PARAMETER_READER_HH_ */ diff --git a/extra_packages/traction-at-split-node-contact/src/common/synchronized_array.hh b/extra_packages/traction-at-split-node-contact/src/common/synchronized_array.hh index 7b8757f76..a88c6ca35 100644 --- a/extra_packages/traction-at-split-node-contact/src/common/synchronized_array.hh +++ b/extra_packages/traction-at-split-node-contact/src/common/synchronized_array.hh @@ -1,204 +1,204 @@ /** * @file synchronized_array.hh * * @author David Simon Kammer * * @date creation: Fri Mar 16 2018 * @date last modification: Tue Sep 29 2020 * * @brief synchronized array: a array can be registered to another (hereafter * called top) array. If an element is added to or removed from the top array, * the registered array removes or adds at the same position an element. The two * arrays stay therefore synchronized. * * * @section LICENSE * * Copyright (©) 2015-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu 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 of the License, or (at your option) any * later version. - * + * * Akantu 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 Lesser General Public License for more * details. - * + * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef AST_SYNCHRONIZED_ARRAY_HH_ #define AST_SYNCHRONIZED_ARRAY_HH_ /* -------------------------------------------------------------------------- */ // std #include // akantu #include "aka_array.hh" namespace akantu { /* -------------------------------------------------------------------------- */ enum SyncChoice { _added, _deleted }; /* -------------------------------------------------------------------------- */ class SynchronizedArrayBase { public: SynchronizedArrayBase() = default; ~SynchronizedArrayBase() = default; virtual ID getID() const { return "call should be virtual"; }; virtual UInt syncDeletedElements(std::vector & delete_el) = 0; virtual UInt syncAddedElements(UInt nb_added_el) = 0; }; /* -------------------------------------------------------------------------- */ template class SynchronizedArray : public SynchronizedArrayBase, protected Array { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: using value_type = typename Array::value_type; using reference = typename Array::reference; using pointer_type = typename Array::pointer_type; using const_reference = typename Array::const_reference; SynchronizedArray(UInt size, UInt nb_component, const_reference value, const ID & id, const_reference default_value, const std::string & restart_name); ~SynchronizedArray() override = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// push_back template inline void push_back(P && value); /// erase inline void erase(UInt i); // template // inline void erase(const iterator & it); /// synchronize elements void syncElements(SyncChoice sync_choice); /// dump restart file void dumpRestartFile(const std::string & file_name) const; /// read restart file void readRestartFile(const std::string & file_name); /// register depending array void registerDependingArray(SynchronizedArrayBase & array); /// function to print the contain of the class void printself(std::ostream & stream, int indent = 0) const override; /// find position of element Int find(const T & elem) const { return Array::find(elem); }; /// set values to zero inline void zero() { Array::zero(); }; // inline void clear() { memset(values, 0, size*nb_component*sizeof(T)); }; /// set all entries of the array to the value t /// @param t value to fill the array with inline void set(T t) { Array::set(t); } /// set template