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;
}
}
}
/* -------------------------------------------------------------------------- */