diff --git a/src/model/phase_field/phasefield.cc b/src/model/phase_field/phasefield.cc index 94ff1f5fd..bd76c3c68 100644 --- a/src/model/phase_field/phasefield.cc +++ b/src/model/phase_field/phasefield.cc @@ -1,409 +1,411 @@ /** * Copyright (©) 2020-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This file is part of Akantu * * 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 "phasefield.hh" #include "aka_common.hh" #include "phase_field_model.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ PhaseField::PhaseField(PhaseFieldModel & model, const ID & id) : Parsable(ParserType::_phasefield, id), id(id), fem(model.getFEEngine()), model(model), g_c("g_c", *this), spatial_dimension(this->model.getSpatialDimension()), element_filter("element_filter", id), damage_on_qpoints("damage", *this), gradd("grad_d", *this), phi("phi", *this), strain("strain", *this), driving_force("driving_force", *this), driving_energy("driving_energy", *this), damage_energy("damage_energy", *this), damage_energy_density("damage_energy_density", *this), dissipated_energy("dissipated_energy", *this) { AKANTU_DEBUG_IN(); /// for each connectivity types allocate the element filer array of the /// material element_filter.initialize(model.getMesh(), _spatial_dimension = spatial_dimension, _element_kind = _ek_regular); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ PhaseField::PhaseField(PhaseFieldModel & model, Int dim, const Mesh & mesh, FEEngine & fe_engine, const ID & id) : Parsable(ParserType::_phasefield, id), id(id), fem(fe_engine), model(model), g_c("g_c", *this), spatial_dimension(this->model.getSpatialDimension()), element_filter("element_filter", id), damage_on_qpoints("damage", *this, dim, fe_engine, this->element_filter), gradd("grad_d", *this, dim, fe_engine, this->element_filter), phi("phi", *this, dim, fe_engine, this->element_filter), strain("strain", *this, dim, fe_engine, this->element_filter), driving_force("driving_force", *this, dim, fe_engine, this->element_filter), driving_energy("driving_energy", *this, dim, fe_engine, this->element_filter), damage_energy("damage_energy", *this, dim, fe_engine, this->element_filter), damage_energy_density("damage_energy_density", *this, dim, fe_engine, this->element_filter), dissipated_energy("dissipated_energy", *this, dim, fe_engine, this->element_filter) { AKANTU_DEBUG_IN(); /// for each connectivity types allocate the element filer array of the /// material element_filter.initialize(mesh, _spatial_dimension = spatial_dimension, _element_kind = _ek_regular); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ PhaseField::~PhaseField() = default; /* -------------------------------------------------------------------------- */ void PhaseField::initialize() { registerParam("name", name, std::string(), _pat_parsable | _pat_readable); registerParam("l0", l0, Real(0.), _pat_parsable | _pat_readable, "length scale parameter"); registerParam("gc", g_c, _pat_parsable | _pat_readable, "critical local fracture energy density"); registerParam("E", E, _pat_parsable | _pat_readable, "Young's modulus"); registerParam("nu", nu, _pat_parsable | _pat_readable, "Poisson ratio"); + registerParam("isotropic", isotropic, true, _pat_parsable | _pat_readable, + "Use isotropic formulation"); damage_on_qpoints.initialize(1); phi.initialize(1); driving_force.initialize(1); driving_energy.initialize(spatial_dimension); gradd.initialize(spatial_dimension); g_c.initialize(1); strain.initialize(spatial_dimension * spatial_dimension); dissipated_energy.initialize(1); damage_energy_density.initialize(1); damage_energy.initialize(spatial_dimension * spatial_dimension); } /* -------------------------------------------------------------------------- */ void PhaseField::initPhaseField() { AKANTU_DEBUG_IN(); this->phi.initializeHistory(); this->resizeInternals(); updateInternalParameters(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseField::resizeInternals() { AKANTU_DEBUG_IN(); for (auto it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) { it->second->resize(); } for (auto it = internal_vectors_int.begin(); it != internal_vectors_int.end(); ++it) { it->second->resize(); } for (auto it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) { it->second->resize(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseField::updateInternalParameters() { this->lambda = this->nu * this->E / ((1 + this->nu) * (1 - 2 * this->nu)); this->mu = this->E / (2 * (1 + this->nu)); } /* -------------------------------------------------------------------------- */ void PhaseField::computeAllDrivingForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); Int spatial_dimension = model.getSpatialDimension(); auto & damage = model.getDamage(); for (const auto & type : element_filter.elementTypes(spatial_dimension, ghost_type)) { auto & elem_filter = element_filter(type, ghost_type); if (elem_filter.empty()) { continue; } // compute the damage on quadrature points auto & damage_interpolated = damage_on_qpoints(type, ghost_type); fem.interpolateOnIntegrationPoints(damage, damage_interpolated, 1, type, ghost_type); auto & gradd_vect = gradd(type, _not_ghost); /// compute @f$\nabla u@f$ fem.gradientOnIntegrationPoints(damage, gradd_vect, 1, type, ghost_type, elem_filter); computeDrivingForce(type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseField::assembleInternalForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); Array & internal_force = model.getInternalForce(); for (auto type : element_filter.elementTypes(_ghost_type = ghost_type)) { auto & elem_filter = element_filter(type, ghost_type); if (elem_filter.empty()) { continue; } auto nb_element = elem_filter.size(); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); auto & driving_force_vect = driving_force(type, ghost_type); Array nt_driving_force(nb_quadrature_points, nb_nodes_per_element); fem.computeNtb(driving_force_vect, nt_driving_force, type, ghost_type, elem_filter); Array int_nt_driving_force(nb_element * nb_quadrature_points, nb_nodes_per_element); fem.integrate(nt_driving_force, int_nt_driving_force, nb_nodes_per_element, type, ghost_type, elem_filter); // damage_energy_on_qpoints = gc*l0 = scalar auto & driving_energy_vect = driving_energy(type, ghost_type); Array bt_driving_energy(nb_element * nb_quadrature_points, nb_nodes_per_element); fem.computeBtD(driving_energy_vect, bt_driving_energy, type, ghost_type, elem_filter); Array int_bt_driving_energy(nb_element, nb_nodes_per_element); fem.integrate(bt_driving_energy, int_bt_driving_energy, nb_nodes_per_element, type, ghost_type, elem_filter); model.getDOFManager().assembleElementalArrayLocalArray( int_nt_driving_force, internal_force, type, ghost_type, -1, elem_filter); model.getDOFManager().assembleElementalArrayLocalArray( int_bt_driving_energy, internal_force, type, ghost_type, -1, elem_filter); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseField::assembleStiffnessMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix"); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { auto & elem_filter = element_filter(type, ghost_type); if (elem_filter.empty()) { AKANTU_DEBUG_OUT(); return; } auto nb_element = elem_filter.size(); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); auto nt_b_n = std::make_unique>( nb_element * nb_quadrature_points, nb_nodes_per_element * nb_nodes_per_element, "N^t*b*N"); auto bt_d_b = std::make_unique>( nb_element * nb_quadrature_points, nb_nodes_per_element * nb_nodes_per_element, "B^t*D*B"); // damage_energy_density_on_qpoints = gc/l0 + phi = scalar auto & damage_energy_density_vect = damage_energy_density(type, ghost_type); // damage_energy_on_qpoints = gc*l0 = scalar auto & damage_energy_vect = damage_energy(type, ghost_type); fem.computeBtDB(damage_energy_vect, *bt_d_b, 2, type, ghost_type, elem_filter); fem.computeNtbN(damage_energy_density_vect, *nt_b_n, type, ghost_type, elem_filter); /// compute @f$ K_{\grad d} = \int_e \mathbf{N}^t * \mathbf{w} * /// \mathbf{N}@f$ auto K_n = std::make_unique>( nb_element, nb_nodes_per_element * nb_nodes_per_element, "K_n"); fem.integrate(*nt_b_n, *K_n, nb_nodes_per_element * nb_nodes_per_element, type, ghost_type, elem_filter); model.getDOFManager().assembleElementalMatricesToMatrix( "K", "damage", *K_n, type, _not_ghost, _symmetric, elem_filter); /// compute @f$ K_{\grad d} = \int_e \mathbf{B}^t * \mathbf{W} * /// \mathbf{B}@f$ auto K_b = std::make_unique>( nb_element, nb_nodes_per_element * nb_nodes_per_element, "K_b"); fem.integrate(*bt_d_b, *K_b, nb_nodes_per_element * nb_nodes_per_element, type, ghost_type, elem_filter); model.getDOFManager().assembleElementalMatricesToMatrix( "K", "damage", *K_b, type, _not_ghost, _symmetric, elem_filter); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseField::computeDissipatedEnergyByElements() { AKANTU_DEBUG_IN(); const Array & damage = model.getDamage(); for (auto type : element_filter.elementTypes(spatial_dimension, _not_ghost)) { Array & elem_filter = element_filter(type, _not_ghost); if (elem_filter.empty()) { continue; } Array & damage_interpolated = damage_on_qpoints(type, _not_ghost); // compute the damage on quadrature points fem.interpolateOnIntegrationPoints(damage, damage_interpolated, 1, type, _not_ghost); Array & gradd_vect = gradd(type, _not_ghost); /// compute @f$\nabla u@f$ fem.gradientOnIntegrationPoints(damage, gradd_vect, 1, type, _not_ghost, elem_filter); computeDissipatedEnergy(type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseField::computeDissipatedEnergy(ElementType /*unused*/) { AKANTU_DEBUG_IN(); AKANTU_TO_IMPLEMENT(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real PhaseField::getEnergy() { AKANTU_DEBUG_IN(); Real edis = 0.; computeDissipatedEnergyByElements(); /// integrate the dissipated energy for each type of elements for (auto type : element_filter.elementTypes(spatial_dimension, _not_ghost)) { edis += fem.integrate(dissipated_energy(type, _not_ghost), type, _not_ghost, element_filter(type, _not_ghost)); } AKANTU_DEBUG_OUT(); return edis; } /* -------------------------------------------------------------------------- */ Real PhaseField::getEnergy(ElementType type, Idx index) { Real edis = 0.; Vector edis_on_quad_points(fem.getNbIntegrationPoints(type)); computeDissipatedEnergyByElement(type, index, edis_on_quad_points); edis = fem.integrate(edis_on_quad_points, Element{type, index, _not_ghost}); return edis; } /* -------------------------------------------------------------------------- */ Real PhaseField::getEnergy(const Element & element) { return getEnergy(element.type, element.element); } /* -------------------------------------------------------------------------- */ void PhaseField::beforeSolveStep() { this->savePreviousState(); this->computeAllDrivingForces(_not_ghost); } /* -------------------------------------------------------------------------- */ void PhaseField::afterSolveStep() {} /* -------------------------------------------------------------------------- */ void PhaseField::savePreviousState() { AKANTU_DEBUG_IN(); for (auto pair : internal_vectors_real) { if (pair.second->hasHistory()) { pair.second->saveCurrentValues(); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseField::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); std::string type = getID().substr(getID().find_last_of(':') + 1); stream << space << "PhaseField Material " << type << " [" << std::endl; Parsable::printself(stream, indent); stream << space << "]" << std::endl; } } // namespace akantu diff --git a/src/model/phase_field/phasefield.hh b/src/model/phase_field/phasefield.hh index 85240799f..95b46d586 100644 --- a/src/model/phase_field/phasefield.hh +++ b/src/model/phase_field/phasefield.hh @@ -1,357 +1,360 @@ /** * Copyright (©) 2020-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This file is part of Akantu * * 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_factory.hh" #include "data_accessor.hh" #include "integration_point.hh" #include "parsable.hh" #include "parser.hh" /* -------------------------------------------------------------------------- */ #include "internal_field.hh" #include "random_internal_field.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_PHASEFIELD_HH__ #define __AKANTU_PHASEFIELD_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class Model; class PhaseFieldModel; class PhaseField; } // namespace akantu namespace akantu { template using InternalPhaseField = InternalFieldTmpl; template <> inline void ParameterTyped>::setAuto( const ParserParameter & in_param) { Parameter::setAuto(in_param); RandomParameter random_param = in_param; param.setRandomDistribution(random_param); } using PhaseFieldFactory = Factory; class PhaseField : public DataAccessor, public Parsable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: PhaseField(const PhaseField & phase) = delete; PhaseField & operator=(const PhaseField & phase) = delete; /// Initialize phasefield with defaults PhaseField(PhaseFieldModel & model, const ID & id = ""); /// Initialize phasefield with custom mesh & fe_engine PhaseField(PhaseFieldModel & model, Int dim, const Mesh & mesh, FEEngine & fe_engine, const ID & id = ""); /// Destructor ~PhaseField() override; protected: void initialize(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: template void registerInternal(InternalPhaseField & /*vect*/) { AKANTU_TO_IMPLEMENT(); } template void unregisterInternal(InternalPhaseField & /*vect*/) { AKANTU_TO_IMPLEMENT(); } /// initialize the phasefield computed parameter virtual void initPhaseField(); /// virtual void beforeSolveStep(); /// virtual void afterSolveStep(); /// assemble the residual for this phasefield virtual void assembleInternalForces(GhostType ghost_type); /// assemble the stiffness matrix for this phasefield virtual void assembleStiffnessMatrix(GhostType ghost_type); /// compute the driving force for this phasefield virtual void computeAllDrivingForces(GhostType ghost_type = _not_ghost); /// save the phi in the phi internal field if needed virtual void savePreviousState(); /// add an element to the local mesh filter inline Int addElement(const Element & element); /// function to print the contain of the class void printself(std::ostream & stream, int indent = 0) const override; protected: /// compute the dissipated energy by element void computeDissipatedEnergyByElements(); /// add an element to the local mesh filter inline Int addElement(const ElementType & type, Idx element, const GhostType & ghost_type); /// resize the internals arrrays virtual void resizeInternals(); /// function called to updatet the internal parameters when the /// modifiable parameters are modified virtual void updateInternalParameters(); // constitutive law for driving force virtual void computeDrivingForce(ElementType /* el_type */, GhostType /* ghost_type */ = _not_ghost) { AKANTU_TO_IMPLEMENT(); } /// compute the dissiapted energy virtual void computeDissipatedEnergy(ElementType el_type); /// compute the dissipated energy for an element virtual void computeDissipatedEnergyByElement(const Element & /*element*/, Vector & /*edis_on_quad_points*/) { AKANTU_TO_IMPLEMENT(); } /// compute the dissipated energy for an element virtual void computeDissipatedEnergyByElement(ElementType /*type*/, Idx /*index*/, Vector & /*edis_on_quad_points*/) { AKANTU_TO_IMPLEMENT(); } /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ public: inline Int 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; template inline void packElementDataHelper(const ElementTypeMapArray & data_to_pack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id = ID()) const; template inline void unpackElementDataHelper(ElementTypeMapArray & data_to_unpack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id = ID()); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ protected: /// return the damage energyfor the provided element virtual Real getEnergy(ElementType type, Idx index); public: /// return the damage energyfor the subset of elements contained /// by the phasefield virtual Real getEnergy(); /// Compute dissipated energy for an individual element Real getEnergy(const Element & element); AKANTU_GET_MACRO(Name, name, const std::string &); AKANTU_GET_MACRO(Model, model, const PhaseFieldModel &) AKANTU_GET_MACRO(ID, id, const ID &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Strain, strain, Real); AKANTU_GET_MACRO(Strain, strain, const ElementTypeMapArray &); AKANTU_GET_MACRO_NOT_CONST(Strain, strain, ElementTypeMapArray &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Damage, damage_on_qpoints, Real); AKANTU_GET_MACRO_NOT_CONST(Damage, damage_on_qpoints, ElementTypeMapArray &); AKANTU_GET_MACRO(Damage, damage_on_qpoints, const ElementTypeMapArray &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementFilter, element_filter, Idx); AKANTU_GET_MACRO(ElementFilter, element_filter, const ElementTypeMapArray &); template const Array & getArray(const ID & id, ElementType type, GhostType ghost_type = _not_ghost) const; template Array & getArray(const ID & id, ElementType type, GhostType ghost_type = _not_ghost); template const InternalPhaseField & getInternal(const ID & id) const; template InternalPhaseField & getInternal(const ID & id); template inline bool isInternal(const ID & id, const ElementKind & element_kind) const; template inline void setParam(const ID & param, T value); inline const Parameter & getParam(const ID & param) const; template void flattenInternal(const std::string & field_id, ElementTypeMapArray & internal_flat, GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined) const; template void inflateInternal(const std::string & field_id, const ElementTypeMapArray & field, GhostType ghost_type = _not_ghost, ElementKind element_kind = _ek_not_defined); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// boolean to know if the material has been initialized bool is_init; std::map *> internal_vectors_real; std::map *> internal_vectors_int; std::map *> internal_vectors_bool; protected: ID id; /// Link to the fem object in the model FEEngine & fem; /// phasefield name std::string name; /// The model to whch the phasefield belong PhaseFieldModel & model; /// length scale parameter Real l0; /// critical energy release rate // Real g_c; RandomInternalField g_c; /// Young's modulus Real E; /// Poisson ratio Real nu; + /// Isotropic formulation + bool isotropic{true}; + /// Lame's first parameter Real lambda; /// Lame's second paramter Real mu; /// spatial dimension Int spatial_dimension; /// list of element handled by the phasefield ElementTypeMapArray element_filter; /// damage arrays ordered by element types InternalPhaseField damage_on_qpoints; /// grad_d arrays ordered by element types InternalPhaseField gradd; /// phi arrays ordered by element types InternalPhaseField phi; /// strain arrays ordered by element types InternalPhaseField strain; /// driving force ordered by element types InternalPhaseField driving_force; /// driving energy ordered by element types InternalPhaseField driving_energy; /// damage energy ordered by element types InternalPhaseField damage_energy; /// damage energy density ordered by element types InternalPhaseField damage_energy_density; /// dissipated energy by element InternalPhaseField dissipated_energy; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const PhaseField & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "phasefield_inline_impl.hh" #include "internal_field_tmpl.hh" #include "random_internal_field_tmpl.hh" /* -------------------------------------------------------------------------- */ #define PHASEFIELD_DEFAULT_ALLOCATOR(id, phase_name) \ [](const ID &, PhaseFieldModel & model, \ const ID & id) -> std::unique_ptr { \ return std::make_unique(model, id); \ } #define INSTANTIATE_PHASEFIELD(id, phase_name) \ static bool phasefield_is_allocated_##id [[gnu::unused]] = \ PhaseFieldFactory::getInstance().registerAllocator( \ #id, PHASEFIELD_DEFAULT_ALLOCATOR(id, phase_name)) #endif diff --git a/src/model/phase_field/phasefields/phasefield_exponential.cc b/src/model/phase_field/phasefields/phasefield_exponential.cc index 3f72ea0de..653ab3e07 100644 --- a/src/model/phase_field/phasefields/phasefield_exponential.cc +++ b/src/model/phase_field/phasefields/phasefield_exponential.cc @@ -1,136 +1,153 @@ /** * Copyright (©) 2020-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This file is part of Akantu * * 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 "phasefield_exponential.hh" #include "aka_common.hh" #include namespace akantu { /* -------------------------------------------------------------------------- */ PhaseFieldExponential::PhaseFieldExponential(PhaseFieldModel & model, const ID & id) : PhaseField(model, id) {} /* -------------------------------------------------------------------------- */ void PhaseFieldExponential::updateInternalParameters() { PhaseField::updateInternalParameters(); for (const auto & type : element_filter.elementTypes(spatial_dimension, _not_ghost)) { for (auto && tuple : zip(make_view(this->damage_energy(type, _not_ghost), spatial_dimension, spatial_dimension), this->g_c(type, _not_ghost))) { Matrix d = Matrix::Identity(spatial_dimension, spatial_dimension) * std::get<1>(tuple) * this->l0; std::get<0>(tuple) = d; } } } /* -------------------------------------------------------------------------- */ void PhaseFieldExponential::computeDrivingForce(ElementType el_type, GhostType ghost_type) { + + if (this->isotropic) { + for (auto && tuple : zip(this->phi(el_type, ghost_type), + this->phi.previous(el_type, ghost_type), + make_view(this->strain(el_type, ghost_type), + spatial_dimension, spatial_dimension))) { + auto & phi_quad = std::get<0>(tuple); + auto & phi_hist_quad = std::get<1>(tuple); + auto & strain = std::get<2>(tuple); + computePhiIsotropicOnQuad(strain, phi_quad, phi_hist_quad); + } + } else { + for (auto && tuple : zip(this->phi(el_type, ghost_type), + this->phi.previous(el_type, ghost_type), + make_view(this->strain(el_type, ghost_type), + spatial_dimension, spatial_dimension))) { + auto & phi_quad = std::get<0>(tuple); + auto & phi_hist_quad = std::get<1>(tuple); + auto & strain = std::get<2>(tuple); + computePhiOnQuad(strain, phi_quad, phi_hist_quad); + } + } + for (auto && tuple : zip(this->phi(el_type, ghost_type), - this->phi.previous(el_type, ghost_type), this->driving_force(el_type, ghost_type), this->damage_energy_density(el_type, ghost_type), - make_view(this->strain(el_type, ghost_type), spatial_dimension, - spatial_dimension), this->damage_on_qpoints(el_type, _not_ghost), make_view(this->driving_energy(el_type, ghost_type), spatial_dimension), make_view(this->damage_energy(el_type, ghost_type), spatial_dimension, spatial_dimension), make_view(this->gradd(el_type, ghost_type), spatial_dimension), this->g_c(el_type, ghost_type))) { - auto & phi_quad = std::get<0>(tuple); - auto & phi_hist_quad = std::get<1>(tuple); - auto & driving_force_quad = std::get<2>(tuple); - auto & dam_energy_density_quad = std::get<3>(tuple); - auto & strain = std::get<4>(tuple); - auto & dam_on_quad = std::get<5>(tuple); - auto & driving_energy_quad = std::get<6>(tuple); - auto & damage_energy_quad = std::get<7>(tuple); - auto & gradd_quad = std::get<8>(tuple); - auto & g_c_quad = std::get<9>(tuple); - - computePhiOnQuad(strain, phi_quad, phi_hist_quad); + auto & driving_force_quad = std::get<1>(tuple); + auto & dam_energy_density_quad = std::get<2>(tuple); + auto & dam_on_quad = std::get<3>(tuple); + auto & driving_energy_quad = std::get<4>(tuple); + auto & damage_energy_quad = std::get<5>(tuple); + auto & gradd_quad = std::get<6>(tuple); + auto & g_c_quad = std::get<7>(tuple); + computeDamageEnergyDensityOnQuad(phi_quad, dam_energy_density_quad, g_c_quad); driving_force_quad = dam_on_quad * dam_energy_density_quad - 2 * phi_quad; driving_energy_quad = damage_energy_quad * gradd_quad; } } + /* -------------------------------------------------------------------------- */ void PhaseFieldExponential::computeDissipatedEnergy(ElementType el_type) { AKANTU_DEBUG_IN(); for (auto && tuple : zip(this->dissipated_energy(el_type, _not_ghost), this->damage_on_qpoints(el_type, _not_ghost), make_view(this->gradd(el_type, _not_ghost), spatial_dimension), this->g_c(el_type, _not_ghost))) { this->computeDissipatedEnergyOnQuad(std::get<1>(tuple), std::get<2>(tuple), std::get<0>(tuple), std::get<3>(tuple)); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseFieldExponential::computeDissipatedEnergyByElement( ElementType type, Idx index, Vector & edis_on_quad_points) { auto gradd_it = this->gradd(type).begin(spatial_dimension); auto gradd_end = this->gradd(type).begin(spatial_dimension); auto damage_it = this->damage_on_qpoints(type).begin(); auto g_c_it = this->g_c(type).begin(); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type); gradd_it += index * nb_quadrature_points; gradd_end += (index + 1) * nb_quadrature_points; damage_it += index * nb_quadrature_points; g_c_it += index * nb_quadrature_points; Real * edis_quad = edis_on_quad_points.data(); for (; gradd_it != gradd_end; ++gradd_it, ++damage_it, ++edis_quad) { this->computeDissipatedEnergyOnQuad(*damage_it, *gradd_it, *edis_quad, *g_c_it); } } void PhaseFieldExponential::computeDissipatedEnergyByElement( const Element & element, Vector & edis_on_quad_points) { computeDissipatedEnergyByElement(element.type, element.element, edis_on_quad_points); } INSTANTIATE_PHASEFIELD(exponential, PhaseFieldExponential); } // namespace akantu diff --git a/src/model/phase_field/phasefields/phasefield_exponential.hh b/src/model/phase_field/phasefields/phasefield_exponential.hh index 3f22bb586..aa7db80f0 100644 --- a/src/model/phase_field/phasefields/phasefield_exponential.hh +++ b/src/model/phase_field/phasefields/phasefield_exponential.hh @@ -1,83 +1,87 @@ /** * Copyright (©) 2020-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This file is part of Akantu * * 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 "phasefield.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_PHASEFIELD_EXPONENTIAL_HH__ #define __AKANTU_PHASEFIELD_EXPONENTIAL_HH__ namespace akantu { class PhaseFieldExponential : public PhaseField { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: PhaseFieldExponential(PhaseFieldModel & model, const ID & id = ""); ~PhaseFieldExponential() override = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// compute the dissiapted energy void computeDissipatedEnergy(ElementType el_type) override; void computeDissipatedEnergyByElement(const Element & element, Vector & edis_on_quad_points) override; protected: void computeDissipatedEnergyByElement(ElementType type, Idx index, Vector & edis_on_quad_points) override; - void computePhiOnQuad(const Matrix & /*strain_quad*/, + inline void computePhiOnQuad(const Matrix & /*strain_quad*/, Real & /*phi_quad*/, Real & /*phi_hist_quad*/); + inline void computePhiIsotropicOnQuad(const Matrix & /*strain_quad*/, + Real & /*phi_quad*/, + Real & /*phi_hist_quad*/); + void computeDrivingForce(ElementType /*el_type*/, GhostType /*ghost_type*/) override; inline void computeDrivingForceOnQuad(const Real & /*phi_quad*/, Real & /*driving_force_quad*/); inline void computeDamageEnergyDensityOnQuad(const Real & /*phi_quad*/, Real & /*dam_energy_quad*/, const Real & /*g_c_quad*/); inline void computeDissipatedEnergyOnQuad(const Real & /*dam_quad*/, const Vector & /*grad_d_quad */, Real & /*energy*/, Real & /*g_c_quad*/); public: void updateInternalParameters() override; }; } // namespace akantu /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "phasefield_exponential_inline_impl.hh" #endif diff --git a/src/model/phase_field/phasefields/phasefield_exponential_inline_impl.hh b/src/model/phase_field/phasefields/phasefield_exponential_inline_impl.hh index b58138ce2..4c8835a75 100644 --- a/src/model/phase_field/phasefields/phasefield_exponential_inline_impl.hh +++ b/src/model/phase_field/phasefields/phasefield_exponential_inline_impl.hh @@ -1,87 +1,101 @@ /** * Copyright (©) 2022-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This file is part of Akantu * * 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 "phasefield_exponential.hh" namespace akantu { inline void PhaseFieldExponential::computeDissipatedEnergyOnQuad( const Real & dam, const Vector & grad_d, Real & edis, Real & g_c_quad) { edis = 0.; for (auto i : arange(spatial_dimension)) { edis += 0.5 * g_c_quad * this->l0 * grad_d[i] * grad_d[i]; } edis += g_c_quad * dam * dam / (2 * this->l0); } /* -------------------------------------------------------------------------- */ inline void PhaseFieldExponential::computeDamageEnergyDensityOnQuad( const Real & phi_quad, Real & dam_energy_quad, const Real & g_c_quad) { dam_energy_quad = 2.0 * phi_quad + g_c_quad / this->l0; } /* -------------------------------------------------------------------------- */ inline void PhaseFieldExponential::computePhiOnQuad(const Matrix & strain_quad, Real & phi_quad, Real & phi_hist_quad) { Matrix strain_plus(spatial_dimension, spatial_dimension); Matrix strain_dir(spatial_dimension, spatial_dimension); Matrix strain_diag_plus(spatial_dimension, spatial_dimension); Vector strain_values(spatial_dimension); Real trace_plus; strain_plus.zero(); strain_dir.zero(); strain_values.zero(); strain_diag_plus.zero(); strain_quad.eig(strain_values, strain_dir); for (Int i = 0; i < spatial_dimension; i++) { strain_diag_plus(i, i) = std::max(Real(0.), strain_values(i)); } Matrix mat_tmp(spatial_dimension, spatial_dimension); Matrix sigma_plus(spatial_dimension, spatial_dimension); mat_tmp = strain_diag_plus * strain_dir.transpose(); strain_plus = strain_dir * mat_tmp; trace_plus = std::max(Real(0.), strain_quad.trace()); for (Int i = 0; i < spatial_dimension; i++) { for (Int j = 0; j < spatial_dimension; j++) { sigma_plus(i, j) = static_cast(i == j) * lambda * trace_plus + 2 * mu * strain_plus(i, j); } } phi_quad = sigma_plus.doubleDot(strain_plus) / 2.; if (phi_quad < phi_hist_quad) { phi_quad = phi_hist_quad; } } +/* -------------------------------------------------------------------------- */ +inline void PhaseFieldExponential::computePhiIsotropicOnQuad( + const Matrix & strain_quad, Real & phi_quad, Real & phi_hist_quad) { + + Real trace = strain_quad.trace(); + + phi_quad = 0.5 * this->lambda * trace * trace + + this->mu * strain_quad.doubleDot(strain_quad); + + if (phi_quad < phi_hist_quad) { + phi_quad = phi_hist_quad; + } +} + } // namespace akantu diff --git a/src/model/solid_mechanics/materials/material_damage/material_phasefield.cc b/src/model/solid_mechanics/materials/material_damage/material_phasefield.cc index 6316b204c..0938837e9 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_phasefield.cc +++ b/src/model/solid_mechanics/materials/material_damage/material_phasefield.cc @@ -1,78 +1,101 @@ /** * Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This file is part of Akantu * * 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_phasefield.hh" #include "aka_common.hh" #include "solid_mechanics_model.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialPhaseField::MaterialPhaseField(SolidMechanicsModel & model, const ID & id) : Parent(model, id), effective_damage("effective_damage", *this) { this->registerParam("eta", eta, Real(0.), _pat_parsable, "eta"); + this->registerParam("is_hybrid", is_hybrid, false, + _pat_parsable | _pat_readable, + "Use hybrid formulation"); this->damage.initialize(0); this->effective_damage.initialize(1); } +/* -------------------------------------------------------------------------- */ template void MaterialPhaseField::computeStress(ElementType el_type, GhostType ghost_type) { - computeEffectiveDamage(el_type, ghost_type); - for (auto && args : getArguments(el_type, ghost_type)) { - computeStressOnQuad(args); + if (this->is_hybrid) { + computeEffectiveDamage(el_type, ghost_type); + + for (auto && args : getArguments(el_type, ghost_type)) { + auto && dam = args["effective_damage"_n]; + computeStressOnQuad(tuple::replace(args, "damage"_n = dam)); + } + } else { + for (auto && args : getArguments(el_type, ghost_type)) { + computeStressOnQuad(args); + } } } /* -------------------------------------------------------------------------- */ template void MaterialPhaseField::computeTangentModuli(ElementType el_type, Array & tangent_matrix, GhostType ghost_type) { computeEffectiveDamage(el_type, ghost_type); - for (auto && args : - getArgumentsTangent(tangent_matrix, el_type, ghost_type)) { - computeTangentModuliOnQuad(args); + if (this->is_hybrid) { + computeEffectiveDamage(el_type, ghost_type); + + for (auto && args : + getArgumentsTangent(tangent_matrix, el_type, ghost_type)) { + auto && dam = args["effective_damage"_n]; + computeTangentModuliOnQuad(tuple::replace(args, "damage"_n = dam)); + } + } else { + for (auto && args : + getArgumentsTangent(tangent_matrix, el_type, ghost_type)) { + computeTangentModuliOnQuad(args); + } } } + /* -------------------------------------------------------------------------- */ template void MaterialPhaseField::computeEffectiveDamage(ElementType el_type, GhostType ghost_type) { for (auto && args : getArguments(el_type, ghost_type)) { computeEffectiveDamageOnQuad(args); } } /* -------------------------------------------------------------------------- */ template class MaterialPhaseField<1>; template class MaterialPhaseField<2>; template class MaterialPhaseField<3>; static bool material_is_allocated_phasefield = instantiateMaterial("phasefield"); } // namespace akantu diff --git a/src/model/solid_mechanics/materials/material_damage/material_phasefield.hh b/src/model/solid_mechanics/materials/material_damage/material_phasefield.hh index 63d55cf1b..769aa451c 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_phasefield.hh +++ b/src/model/solid_mechanics/materials/material_damage/material_phasefield.hh @@ -1,101 +1,105 @@ /** * Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This file is part of Akantu * * 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" #include "material_damage.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_PHASEFIELD_HH__ #define __AKANTU_MATERIAL_PHASEFIELD_HH__ namespace akantu { template class MaterialPhaseField : public MaterialDamage { using Parent = MaterialDamage; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialPhaseField(SolidMechanicsModel & model, const ID & id = ""); ~MaterialPhaseField() override = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// constitutive law for all element of a type void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost) override; /// compute the tangent stiffness matrix for an element type void computeTangentModuli(ElementType el_type, Array & tangent_matrix, GhostType ghost_type = _not_ghost) override; /* ------------------------------------------------------------------------ */ decltype(auto) getArguments(ElementType el_type, GhostType ghost_type = _not_ghost) { return zip_append(Parent::getArguments(el_type, ghost_type), "effective_damage"_n = make_view( this->effective_damage(el_type, ghost_type))); } decltype(auto) getArgumentsTangent(Array & tangent_matrix, ElementType el_type, GhostType ghost_type) { return zip_append( Parent::getArgumentsTangent(tangent_matrix, el_type, ghost_type), "effective_damage"_n = make_view(this->effective_damage(el_type, ghost_type))); } protected: /// constitutive law for a given quadrature point template inline void computeStressOnQuad(Args && args); /// compute the tangent stiffness matrix for a given quadrature point template inline void computeTangentModuliOnQuad(Args && args); /// Compute the effective damage void computeEffectiveDamage(ElementType el_type, GhostType ghost_type = _not_ghost); template inline void computeEffectiveDamageOnQuad(Args && args); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: + /// Residual stiffness parameter Real eta; + /// Phasefield isotropic + bool is_hybrid; + // effective damage to conserve stiffness in compression InternalField effective_damage; }; } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_phasefield_inline_impl.hh" #endif /* __AKANTU_MATERIAL_PHASEFIELD_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_damage/material_phasefield_inline_impl.hh b/src/model/solid_mechanics/materials/material_damage/material_phasefield_inline_impl.hh index 2db8c6240..f732948fa 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_phasefield_inline_impl.hh +++ b/src/model/solid_mechanics/materials/material_damage/material_phasefield_inline_impl.hh @@ -1,89 +1,89 @@ /** * Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This file is part of Akantu * * 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_phasefield.hh" #ifndef AKANTU_MATERIAL_PHASEFIELD_INLINE_IMPL_HH_ #define AKANTU_MATERIAL_PHASEFIELD_INLINE_IMPL_HH_ /* -------------------------------------------------------------------------- */ namespace akantu { template template inline void MaterialPhaseField::computeStressOnQuad(Args && args) { MaterialElastic::computeStressOnQuad(args); - auto && dam = args["effective_damage"_n]; + auto && dam = args["damage"_n]; args["sigma"_n] *= (1 - dam) * (1 - dam) + eta; } /* -------------------------------------------------------------------------- */ template template void MaterialPhaseField::computeTangentModuliOnQuad(Args && args) { MaterialElastic::computeTangentModuliOnQuad(args); - auto dam = args["effective_damage"_n]; + auto dam = args["damage"_n]; args["tangent_moduli"_n] *= (1 - dam) * (1 - dam) + eta; } /* -------------------------------------------------------------------------- */ template template inline void MaterialPhaseField::computeEffectiveDamageOnQuad(Args && args) { using Mat = Matrix; auto strain = Material::gradUToEpsilon(args["grad_u"_n]); Mat strain_dir; Vector strain_values; strain.eig(strain_values, strain_dir); Mat strain_diag_plus; Mat strain_diag_minus; strain_diag_plus.zero(); strain_diag_minus.zero(); for (UInt i = 0; i < dim; i++) { strain_diag_plus(i, i) = std::max(Real(0.), strain_values(i)); strain_diag_minus(i, i) = std::min(Real(0.), strain_values(i)); } Mat strain_plus = strain_dir * strain_diag_plus * strain_dir.transpose(); Mat strain_minus = strain_dir * strain_diag_minus * strain_dir.transpose(); auto trace_plus = std::max(Real(0.), strain.trace()); auto trace_minus = std::min(Real(0.), strain.trace()); Mat sigma_plus = Mat::Identity() * trace_plus * this->lambda + 2. * this->mu * strain_plus; Mat sigma_minus = Mat::Identity() * trace_minus * this->lambda + 2. * this->mu * strain_minus; auto strain_energy_plus = sigma_plus.doubleDot(strain_plus) / 2.; auto strain_energy_minus = sigma_minus.doubleDot(strain_minus) / 2.; args["effective_damage"_n] = args["damage"_n] * (strain_energy_minus < strain_energy_plus); } } // namespace akantu #endif diff --git a/test/test_model/test_phase_field_model/CMakeLists.txt b/test/test_model/test_phase_field_model/CMakeLists.txt index ac8b16ff0..5386b21be 100644 --- a/test/test_model/test_phase_field_model/CMakeLists.txt +++ b/test/test_model/test_phase_field_model/CMakeLists.txt @@ -1,52 +1,52 @@ #=============================================================================== # Copyright (©) 2010-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # This file is part of Akantu # # 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 . # #=============================================================================== register_test(test_phasefield_selector SOURCES test_phasefield_selector.cc FILES_TO_COPY phasefield_selector.dat phasefield_selector.msh PACKAGE phase_field ) register_test(test_phase_solid_coupling SOURCES test_phase_solid_coupling.cc FILES_TO_COPY material_coupling.dat test_one_element.msh PACKAGE phase_field ) register_test(test_phase_field_anisotropic SOURCES test_phase_field_anisotropic.cc - FILES_TO_COPY material_coupling.dat test_one_element.msh + FILES_TO_COPY material_hybrid.dat test_one_element.msh PACKAGE phase_field ) register_test(test_phase_solid_explicit SOURCES test_phase_solid_explicit.cc FILES_TO_COPY material_coupling.dat test_one_element.msh PACKAGE phase_field ) register_test(test_multi_material SOURCES test_multi_material.cc FILES_TO_COPY material_multiple.dat test_two_element.msh PACKAGE phase_field UNSTABLE ) diff --git a/test/test_model/test_phase_field_model/material_hybrid.dat b/test/test_model/test_phase_field_model/material_hybrid.dat new file mode 100644 index 000000000..cede2f624 --- /dev/null +++ b/test/test_model/test_phase_field_model/material_hybrid.dat @@ -0,0 +1,21 @@ +model solid_mechanics_model [ + material phasefield [ + name = plate + rho = 1. + E = 210.0 + nu = 0.3 + Plane_Stress = false + is_hybrid = true + ] +] + +model phase_field_model [ + phasefield exponential [ + name = plate + E = 210.0 + nu = 0.3 + gc = 5e-3 + l0 = 0.1 + isotropic = false + ] +] diff --git a/test/test_model/test_phase_field_model/test_phase_field_anisotropic.cc b/test/test_model/test_phase_field_model/test_phase_field_anisotropic.cc index 5f1731554..34c26e976 100644 --- a/test/test_model/test_phase_field_model/test_phase_field_anisotropic.cc +++ b/test/test_model/test_phase_field_model/test_phase_field_anisotropic.cc @@ -1,167 +1,168 @@ /** * Copyright (©) 2021-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * This file is part of Akantu * * 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 "coupler_solid_phasefield.hh" #include "material.hh" #include "material_phasefield.hh" #include "non_linear_solver.hh" #include "phase_field_model.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; const UInt spatial_dimension = 2; /* -------------------------------------------------------------------------- */ void applyDisplacement(SolidMechanicsModel &, Real &); /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { std::ofstream os("data.csv"); os << "#strain stress damage analytical_sigma analytical_damage" << std::endl; - initialize("material_coupling.dat", argc, argv); + initialize("material_hybrid.dat", argc, argv); Mesh mesh(spatial_dimension); mesh.read("test_one_element.msh"); CouplerSolidPhaseField coupler(mesh); auto & model = coupler.getSolidMechanicsModel(); auto & phase = coupler.getPhaseFieldModel(); model.initFull(_analysis_method = _static); auto && selector = std::make_shared>( "physical_names", phase); phase.setPhaseFieldSelector(selector); phase.initFull(_analysis_method = _static); model.setBaseName("phase_solid"); model.addDumpField("stress"); model.addDumpField("grad_u"); model.addDumpFieldVector("displacement"); model.addDumpField("damage"); model.dump(); UInt nbSteps = 1500; Real increment = 1e-4; auto & stress = model.getMaterial(0).getArray("stress", _quadrangle_4); auto & damage = model.getMaterial(0).getArray("damage", _quadrangle_4); Real analytical_damage{0.}; Real analytical_sigma{0.}; auto & phasefield = phase.getPhaseField(0); const Real E = phasefield.getParam("E"); const Real nu = phasefield.getParam("nu"); Real c22 = E * (1 - nu) / ((1 + nu) * (1 - 2 * nu)); const Real gc = phasefield.getParam("gc"); const Real l0 = phasefield.getParam("l0"); Real error_stress{0.}; Real error_damage{0.}; Real max_strain{0.}; for (UInt s = 0; s < nbSteps; ++s) { Real axial_strain{0.}; if (s < 500) { axial_strain = increment * s; } else if (s < 1000) { axial_strain = (1500 - 2 * double(s)) * increment; } else { axial_strain = (3 * double(s) - 3500) * increment; } applyDisplacement(model, axial_strain); if (axial_strain > max_strain) { max_strain = axial_strain; } coupler.solve("static", "static"); analytical_damage = max_strain * max_strain * c22 / (gc / l0 + max_strain * max_strain * c22); if (axial_strain < 0.) { analytical_sigma = c22 * axial_strain; } else { analytical_sigma = c22 * axial_strain * (1 - analytical_damage) * (1 - analytical_damage); } - error_stress = std::abs(analytical_sigma - stress(0, 3)) / analytical_sigma; + error_stress = + std::abs(analytical_sigma - stress(0, 3)) / std::abs(analytical_sigma); error_damage = std::abs(analytical_damage - damage(0)) / analytical_damage; if ((error_damage > 1e-8 or error_stress > 1e-8) and - std::abs(axial_strain) < 1e-13) { + std::abs(axial_strain) > 1e-13) { std::cerr << std::left << std::setw(15) << "Error damage: " << error_damage << std::endl; std::cerr << std::left << std::setw(15) << "Error stress: " << error_stress << std::endl; return EXIT_FAILURE; } os << axial_strain << " " << stress(0, 3) << " " << damage(0) << " " << analytical_sigma << " " << analytical_damage << " " << error_stress << " " << error_damage << std::endl; model.dump(); } os.close(); finalize(); return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */ void applyDisplacement(SolidMechanicsModel & model, Real & increment) { auto & displacement = model.getDisplacement(); auto & positions = model.getMesh().getNodes(); auto & blocked_dofs = model.getBlockedDOFs(); for (Idx n = 0; n < model.getMesh().getNbNodes(); ++n) { if (positions(n, 1) == -0.5) { displacement(n, 0) = 0; displacement(n, 1) = 0; blocked_dofs(n, 0) = true; blocked_dofs(n, 1) = true; } else { displacement(n, 0) = 0; displacement(n, 1) = increment; blocked_dofs(n, 0) = true; blocked_dofs(n, 1) = true; } } } /* -------------------------------------------------------------------------- */