diff --git a/examples/c++/phase_field_model/phase_field_notch.cc b/examples/c++/phase_field_model/phase_field_notch.cc
index d787163c2..e3aead0f9 100644
--- a/examples/c++/phase_field_model/phase_field_notch.cc
+++ b/examples/c++/phase_field_model/phase_field_notch.cc
@@ -1,124 +1,124 @@
/**
* Copyright (©) 2018-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 <http://www.gnu.org/licenses/>.
*/
/* -------------------------------------------------------------------------- */
#include "coupler_solid_phasefield.hh"
#include "phase_field_element_filter.hh"
/* -------------------------------------------------------------------------- */
#include <chrono>
#include <iostream>
/* -------------------------------------------------------------------------- */
using namespace akantu;
/* -------------------------------------------------------------------------- */
using clk = std::chrono::high_resolution_clock;
using second = std::chrono::duration<double>;
using millisecond = std::chrono::duration<double, std::milli>;
const Int spatial_dimension = 2;
/* -------------------------------------------------------------------------- */
int main(int argc, char * argv[]) {
initialize("material_notch.dat", argc, argv);
// create mesh
Mesh mesh(spatial_dimension);
mesh.read("square_notch.msh");
// The model coupler contains the solid mechanics and the phase-field models
CouplerSolidPhaseField coupler(mesh);
auto & model = coupler.getSolidMechanicsModel();
auto & phase = coupler.getPhaseFieldModel();
// Each model can bet set separately
model.initFull(_analysis_method = _static);
phase.initFull(_analysis_method = _static);
// Dirichlet BC
model.applyBC(BC::Dirichlet::FixedValue(0., _y), "bottom");
model.applyBC(BC::Dirichlet::FixedValue(0., _x), "left");
// Dumper settings
model.setBaseName("phase_notch");
model.addDumpField("stress");
model.addDumpField("grad_u");
model.addDumpFieldVector("displacement");
model.addDumpField("damage");
model.dump();
const Int nb_steps = 1000;
- Real increment = 6e-6;
- Int nb_staggered_steps = 5;
+ Real increment = 4e-6;
+ Int nb_staggered_steps = 10;
auto start_time = clk::now();
// Main loop over the loading steps
for (Int s = 1; s < nb_steps; ++s) {
- if (s >= 500) {
- increment = 2e-6;
- nb_staggered_steps = 10;
- }
+ // if (s >= 500) {
+ // increment = 2e-6;
+ // nb_staggered_steps = 10;
+ // }
if (s % 200 == 0) {
constexpr std::array<char, 5> wheel{"/-\\|"};
auto elapsed = clk::now() - start_time;
auto time_per_step = elapsed / s;
int idx = (s / 10) % 4;
std::cout << "\r[" << wheel.at(idx) << "] " << std::setw(5) << s << "/"
<< nb_steps << " (" << std::setprecision(2) << std::fixed
<< std::setw(8) << millisecond(time_per_step).count()
<< "ms/step - elapsed: " << std::setw(8)
<< second(elapsed).count() << "s - ETA: " << std::setw(8)
<< second((nb_steps - s) * time_per_step).count() << "s)"
<< std::string(' ', 20) << std::flush;
}
model.applyBC(BC::Dirichlet::IncrementValue(increment, _y), "top");
/* At each step, the two solvers are called alternately. Here a fixed number
* of staggered iterations is set for simplicity but a convergence based on
* the displacements, damage and/or energy can also be used to exit the
* internal loop.*/
for (Idx i = 0; i < nb_staggered_steps; ++i) {
coupler.solve();
}
if (s % 100 == 0) {
model.dump();
}
}
/* Here a damage threshold is set and a mesh clustering function is called
* using the phase-field element filter. This allows to cluster the mesh into
* groups of elements separated by the mesh boundaries and "broken" elements
* with damage above the threshold. */
Real damage_limit = 0.15;
auto global_nb_clusters = mesh.createClusters(
spatial_dimension, "crack", PhaseFieldElementFilter(phase, damage_limit));
model.dumpGroup("crack_0");
model.dumpGroup("crack_1");
std::cout << std::endl;
std::cout << "Nb clusters: " << global_nb_clusters << std::endl;
finalize();
return EXIT_SUCCESS;
}
diff --git a/examples/c++/phase_field_model/phase_field_parallel.cc b/examples/c++/phase_field_model/phase_field_parallel.cc
index 8f419db65..b697949dd 100644
--- a/examples/c++/phase_field_model/phase_field_parallel.cc
+++ b/examples/c++/phase_field_model/phase_field_parallel.cc
@@ -1,129 +1,132 @@
/**
* @file phase_field_parallel.cc
*
* @author Mohit Pundir <mohit.pundir@ethz.ch>
*
* @date creation: Mon May 09 2022
* @date last modification: Mon May 09 2022
*
* @brief Example of phase field model in parallel
*
*
* @section LICENSE
*
* Copyright (©) 2018-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "coupler_solid_phasefield.hh"
/* -------------------------------------------------------------------------- */
#include <chrono>
#include <fstream>
#include <iostream>
/* -------------------------------------------------------------------------- */
using namespace akantu;
using clk = std::chrono::high_resolution_clock;
using second = std::chrono::duration<double>;
using millisecond = std::chrono::duration<double, std::milli>;
const UInt spatial_dimension = 2;
int main(int argc, char * argv[]) {
initialize("material_notch.dat", argc, argv);
// create mesh
Mesh mesh(spatial_dimension);
const auto & comm = Communicator::getStaticCommunicator();
Int prank = comm.whoAmI();
if (prank == 0) {
// Read the mesh
mesh.read("square_notch.msh");
}
mesh.distribute();
CouplerSolidPhaseField coupler(mesh);
auto & model = coupler.getSolidMechanicsModel();
auto & phase = coupler.getPhaseFieldModel();
model.initFull(_analysis_method = _static);
auto && mat_selector =
std::make_shared<MeshDataMaterialSelector<std::string>>("physical_names",
model);
model.setMaterialSelector(mat_selector);
auto && selector = std::make_shared<MeshDataPhaseFieldSelector<std::string>>(
"physical_names", phase);
phase.setPhaseFieldSelector(selector);
phase.initFull(_analysis_method = _static);
model.applyBC(BC::Dirichlet::FixedValue(0., _y), "bottom");
model.applyBC(BC::Dirichlet::FixedValue(0., _x), "left");
model.setBaseName("phase_notch_parallel");
model.addDumpField("stress");
model.addDumpField("grad_u");
model.addDumpFieldVector("displacement");
model.addDumpField("damage");
if (mesh.isDistributed()) {
- // phase.addDumpField("partitions");
+ phase.addDumpField("partitions");
}
phase.dump();
UInt nbSteps = 1000;
- Real increment = 6e-6;
- UInt nb_staggered_steps = 5;
+ Real increment = 4e-6;
+ UInt nb_staggered_steps = 10;
+
+ Int nbNodes = mesh.getNbNodes();
+ std::cout << "Number of nodes: " << nbNodes << std::endl;
auto start_time = clk::now();
for (UInt s = 1; s < nbSteps; ++s) {
- if (s >= 500) {
- increment = 2e-6;
- nb_staggered_steps = 10;
- }
+ // if (s >= 500) {
+ // increment = 2e-6;
+ // nb_staggered_steps = 10;
+ // }
if (s % 10 == 0 && prank == 0) {
constexpr char wheel[] = "/-\\|";
auto elapsed = clk::now() - start_time;
auto time_per_step = elapsed / s;
std::cout << "\r[" << wheel[(s / 10) % 4] << "] " << std::setw(5) << s
<< "/" << nbSteps << " (" << std::setprecision(2) << std::fixed
<< std::setw(8) << millisecond(time_per_step).count()
<< "ms/step - elapsed: " << std::setw(8)
<< second(elapsed).count() << "s - ETA: " << std::setw(8)
<< second((nbSteps - s) * time_per_step).count() << "s)"
<< std::string(' ', 20) << std::flush;
}
model.applyBC(BC::Dirichlet::IncrementValue(increment, _y), "top");
for (UInt i = 0; i < nb_staggered_steps; ++i) {
coupler.solve();
}
if (s % 100 == 0) {
model.dump();
}
}
return 0;
}
diff --git a/src/model/phase_field/phase_field_model.cc b/src/model/phase_field/phase_field_model.cc
index ba5158b27..9265fff34 100644
--- a/src/model/phase_field/phase_field_model.cc
+++ b/src/model/phase_field/phase_field_model.cc
@@ -1,504 +1,504 @@
/**
* Copyright (©) 2018-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 <http://www.gnu.org/licenses/>.
*/
/* -------------------------------------------------------------------------- */
#include "phase_field_model.hh"
#include "dumpable_inline_impl.hh"
#include "element_synchronizer.hh"
#include "fe_engine_template.hh"
#include "group_manager_inline_impl.hh"
#include "integrator_gauss.hh"
#include "shape_lagrange.hh"
/* -------------------------------------------------------------------------- */
#include "dumper_element_partition.hh"
#include "dumper_elemental_field.hh"
#include "dumper_internal_material_field.hh"
#include "dumper_iohelper_paraview.hh"
#include <algorithm>
#include <utility>
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
PhaseFieldModel::PhaseFieldModel(Mesh & mesh, Int dim, const ID & id,
std::shared_ptr<DOFManager> dof_manager,
ModelType model_type)
: Parent(mesh, model_type, dim, id) {
AKANTU_DEBUG_IN();
this->initDOFManager(std::move(dof_manager));
this->registerFEEngineObject<FEEngineType>("PhaseFieldFEEngine", mesh,
Model::spatial_dimension);
this->mesh.registerDumper<DumperParaview>("phase_field", id, true);
this->mesh.addDumpMesh(mesh, Model::spatial_dimension, _not_ghost,
_ek_regular);
if (this->mesh.isDistributed()) {
auto & synchronizer = this->mesh.getElementSynchronizer();
this->registerSynchronizer(synchronizer, SynchronizationTag::_pfm_damage);
this->registerSynchronizer(synchronizer, SynchronizationTag::_for_dump);
}
this->parser_type = ParserType::_phasefield;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
MatrixType PhaseFieldModel::getMatrixType(const ID & matrix_id) const {
if (matrix_id == "K") {
return _symmetric;
}
return _mt_not_defined;
}
/* -------------------------------------------------------------------------- */
void PhaseFieldModel::initFullImpl(const ModelOptions & options) {
Parent::initFullImpl(options);
this->initBC(*this, *damage, *external_force);
}
/* -------------------------------------------------------------------------- */
void PhaseFieldModel::assembleMatrix(const ID & matrix_id) {
if (matrix_id == "K") {
this->assembleStiffnessMatrix();
} else {
AKANTU_ERROR("Unknown Matrix ID for PhaseFieldModel : " << matrix_id);
}
}
/* -------------------------------------------------------------------------- */
void PhaseFieldModel::predictor() {
// AKANTU_TO_IMPLEMENT();
}
/* -------------------------------------------------------------------------- */
void PhaseFieldModel::corrector() {
// AKANTU_TO_IMPLEMENT();
}
/* -------------------------------------------------------------------------- */
void PhaseFieldModel::initSolver(TimeStepSolverType time_step_solver_type,
NonLinearSolverType /*unused*/) {
DOFManager & dof_manager = this->getDOFManager();
this->allocNodalField(this->damage, 1, "damage");
this->allocNodalField(this->external_force, 1, "external_force");
this->allocNodalField(this->internal_force, 1, "internal_force");
this->allocNodalField(this->blocked_dofs, 1, "blocked_dofs");
this->allocNodalField(this->previous_damage, 1, "previous_damage");
if (!dof_manager.hasDOFs("damage")) {
dof_manager.registerDOFs("damage", *this->damage, _dst_nodal);
dof_manager.registerBlockedDOFs("damage", *this->blocked_dofs);
dof_manager.registerDOFsPrevious("damage", *this->previous_damage);
}
if (time_step_solver_type == TimeStepSolverType::_dynamic) {
AKANTU_TO_IMPLEMENT();
}
}
/* -------------------------------------------------------------------------- */
FEEngine & PhaseFieldModel::getFEEngineBoundary(const ID & name) {
return dynamic_cast<FEEngine &>(getFEEngineClassBoundary<FEEngineType>(name));
}
/* -------------------------------------------------------------------------- */
TimeStepSolverType PhaseFieldModel::getDefaultSolverType() const {
return TimeStepSolverType::_static;
}
/* -------------------------------------------------------------------------- */
std::tuple<ID, TimeStepSolverType>
PhaseFieldModel::getDefaultSolverID(const AnalysisMethod & method) {
switch (method) {
case _explicit_lumped_mass: {
return std::make_tuple("explicit_lumped",
TimeStepSolverType::_dynamic_lumped);
}
case _explicit_consistent_mass: {
return std::make_tuple("explicit", TimeStepSolverType::_dynamic);
}
case _static: {
return std::make_tuple("static", TimeStepSolverType::_static);
}
case _implicit_dynamic: {
return std::make_tuple("implicit", TimeStepSolverType::_dynamic);
}
default:
return std::make_tuple("unknown", TimeStepSolverType::_not_defined);
}
}
/* -------------------------------------------------------------------------- */
ModelSolverOptions PhaseFieldModel::getDefaultSolverOptions(
const TimeStepSolverType & type) const {
ModelSolverOptions options;
switch (type) {
case TimeStepSolverType::_dynamic_lumped: {
options.non_linear_solver_type = NonLinearSolverType::_lumped;
options.integration_scheme_type["damage"] =
IntegrationSchemeType::_central_difference;
options.solution_type["damage"] = IntegrationScheme::_acceleration;
break;
}
case TimeStepSolverType::_static: {
options.non_linear_solver_type = NonLinearSolverType::_newton_raphson;
options.integration_scheme_type["damage"] =
IntegrationSchemeType::_pseudo_time;
options.solution_type["damage"] = IntegrationScheme::_not_defined;
break;
}
case TimeStepSolverType::_dynamic: {
options.non_linear_solver_type = NonLinearSolverType::_newton_raphson;
options.integration_scheme_type["damage"] =
IntegrationSchemeType::_backward_euler;
options.solution_type["damage"] = IntegrationScheme::_damage;
break;
}
default:
AKANTU_EXCEPTION(type << " is not a valid time step solver type");
}
return options;
}
/* -------------------------------------------------------------------------- */
Real PhaseFieldModel::getEnergy(const ID & energy_id) {
AKANTU_DEBUG_IN();
Real energy = 0.;
for_each_constitutive_law([&energy, &energy_id](auto && phase_field) {
energy += phase_field.getEnergy(energy_id);
});
/// reduction sum over all processors
mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum);
AKANTU_DEBUG_OUT();
return energy;
}
/* -------------------------------------------------------------------------- */
Real PhaseFieldModel::getEnergy(const ID & energy_id, const Element & element) {
auto pf_element = element;
auto phase_index = this->getConstitutiveLawByElement()(element);
pf_element.element = this->getConstitutiveLawLocalNumbering()(element);
Real energy =
this->getConstitutiveLaw(phase_index).getEnergy(energy_id, pf_element);
return energy;
}
/* -------------------------------------------------------------------------- */
Real PhaseFieldModel::getEnergy(const ID & energy_id, const ID & group_id) {
auto && group = mesh.getElementGroup(group_id);
auto energy = 0.;
for (auto && type : group.elementTypes()) {
for (auto el : group.getElementsIterable(type)) {
energy += this->getEnergy(energy_id, el);
}
}
/// reduction sum over all processors
mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum);
return energy;
}
/* -------------------------------------------------------------------------- */
void PhaseFieldModel::beforeSolveStep() {
for_each_constitutive_law(
[](auto && phasefield) { phasefield.beforeSolveStep(); });
}
/* -------------------------------------------------------------------------- */
void PhaseFieldModel::afterSolveStep(bool converged) {
if (not converged) {
return;
}
- for (auto && values : zip(*damage, *previous_damage)) {
- auto & dam = std::get<0>(values);
- auto & prev_dam = std::get<1>(values);
+ // for (auto && values : zip(*damage, *previous_damage)) {
+ // auto & dam = std::get<0>(values);
+ // auto & prev_dam = std::get<1>(values);
- prev_dam = dam;
- }
+ // prev_dam = dam;
+ // }
for (auto && values : zip(*damage, *blocked_dofs)) {
auto & dam = std::get<0>(values);
auto & blocked = std::get<1>(values);
dam = std::min(1., dam);
if (!blocked) {
blocked = Math::are_float_equal(dam, 1.);
}
}
for_each_constitutive_law(
[](auto && phasefield) { phasefield.afterSolveStep(); });
}
/* -------------------------------------------------------------------------- */
void PhaseFieldModel::assembleStiffnessMatrix() {
AKANTU_DEBUG_INFO("Assemble the new stiffness matrix");
if (!this->getDOFManager().hasMatrix("K")) {
this->getDOFManager().getNewMatrix("K", getMatrixType("K"));
}
this->getDOFManager().zeroMatrix("K");
for_each_constitutive_law([](auto && phasefield) {
phasefield.assembleStiffnessMatrix(_not_ghost);
});
}
/* -------------------------------------------------------------------------- */
void PhaseFieldModel::assembleResidual() {
this->assembleInternalForces();
this->getDOFManager().assembleToResidual("damage", *this->external_force, 1);
this->getDOFManager().assembleToResidual("damage", *this->internal_force, 1);
}
/* -------------------------------------------------------------------------- */
void PhaseFieldModel::assembleInternalForces() {
AKANTU_DEBUG_INFO("Assemble the internal forces");
this->internal_force->zero();
for_each_constitutive_law([](auto && phasefield) {
phasefield.computeAllDrivingForces(_not_ghost);
});
this->asynchronousSynchronize(SynchronizationTag::_pfm_damage);
// assemble the forces due to local driving forces
AKANTU_DEBUG_INFO("Assemble residual for local elements");
for_each_constitutive_law([](auto && phasefield) {
phasefield.assembleInternalForces(_not_ghost);
});
this->waitEndSynchronize(SynchronizationTag::_pfm_damage);
// assemble the forces due to local driving forces
AKANTU_DEBUG_INFO("Assemble residual for ghost elements");
for_each_constitutive_law(
[](auto && phasefield) { phasefield.assembleInternalForces(_ghost); });
}
/* -------------------------------------------------------------------------- */
void PhaseFieldModel::savePreviousState() {
for_each_constitutive_law(
[](auto && phasefield) { phasefield.savePreviousState(); });
}
/* -------------------------------------------------------------------------- */
void PhaseFieldModel::assembleLumpedMatrix(const ID & /*matrix_id*/) {}
/* -------------------------------------------------------------------------- */
void PhaseFieldModel::setTimeStep(Real time_step, const ID & solver_id) {
Model::setTimeStep(time_step, solver_id);
this->mesh.getDumper("phase_field").setTimeStep(time_step);
}
/* -------------------------------------------------------------------------- */
Int PhaseFieldModel::getNbData(const Array<Element> & elements,
const SynchronizationTag & tag) const {
Int size = 0;
Int nb_nodes_per_element = 0;
for (const Element & el : elements) {
nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type);
}
switch (tag) {
case SynchronizationTag::_for_dump: {
// damage
size += nb_nodes_per_element * sizeof(Real);
break;
}
// case SynchronizationTag::_pfm_damage: {
// size += nb_nodes_per_element * sizeof(Real);
// break;
// }
default: {
// AKANTU_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
size += Parent::getNbData(elements, tag);
return size;
}
/* -------------------------------------------------------------------------- */
void PhaseFieldModel::packData(CommunicationBuffer & buffer,
const Array<Element> & elements,
const SynchronizationTag & tag) const {
switch (tag) {
case SynchronizationTag::_for_dump: {
packNodalDataHelper(*damage, buffer, elements, mesh);
break;
}
// case SynchronizationTag::_pfm_damage: {
// packNodalDataHelper(*damage, buffer, elements, mesh);
// break;
// }
default: {
// AKANTU_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
Parent::packData(buffer, elements, tag);
}
/* -------------------------------------------------------------------------- */
void PhaseFieldModel::unpackData(CommunicationBuffer & buffer,
const Array<Element> & elements,
const SynchronizationTag & tag) {
switch (tag) {
case SynchronizationTag::_for_dump: {
unpackNodalDataHelper(*damage, buffer, elements, mesh);
break;
}
// case SynchronizationTag::_pfm_damage: {
// unpackNodalDataHelper(*damage, buffer, elements, mesh);
// break;
// }
default: {
// AKANTU_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
Parent::unpackData(buffer, elements, tag);
}
/* -------------------------------------------------------------------------- */
Int PhaseFieldModel::getNbData(const Array<Idx> & indexes,
const SynchronizationTag & tag) const {
Int size = 0;
Int nb_nodes = indexes.size();
switch (tag) {
case SynchronizationTag::_for_dump: {
size += nb_nodes * sizeof(Real);
break;
}
case SynchronizationTag::_pfm_damage: {
size += nb_nodes * sizeof(Real);
break;
}
default: {
AKANTU_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
return size;
}
/* -------------------------------------------------------------------------- */
void PhaseFieldModel::packData(CommunicationBuffer & buffer,
const Array<Idx> & indexes,
const SynchronizationTag & tag) const {
switch (tag) {
case SynchronizationTag::_for_dump: {
packDOFDataHelper(*damage, buffer, indexes);
break;
}
case SynchronizationTag::_pfm_damage: {
packDOFDataHelper(*damage, buffer, indexes);
break;
}
default: {
AKANTU_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
}
/* -------------------------------------------------------------------------- */
void PhaseFieldModel::unpackData(CommunicationBuffer & buffer,
const Array<Idx> & indexes,
const SynchronizationTag & tag) {
switch (tag) {
case SynchronizationTag::_for_dump: {
unpackDOFDataHelper(*damage, buffer, indexes);
break;
}
case SynchronizationTag::_pfm_damage: {
unpackDOFDataHelper(*damage, buffer, indexes);
break;
}
default: {
AKANTU_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
}
/* -------------------------------------------------------------------------- */
std::shared_ptr<dumpers::Field>
PhaseFieldModel::createNodalFieldBool(const std::string & field_name,
const std::string & group_name,
bool /*unused*/) {
std::map<std::string, Array<bool> *> uint_nodal_fields;
uint_nodal_fields["blocked_dofs"] = blocked_dofs.get();
return mesh.createNodalField(uint_nodal_fields[field_name], group_name);
std::shared_ptr<dumpers::Field> field;
return field;
}
/* -------------------------------------------------------------------------- */
std::shared_ptr<dumpers::Field>
PhaseFieldModel::createNodalFieldReal(const std::string & field_name,
const std::string & group_name,
bool /*unused*/) {
std::map<std::string, Array<Real> *> real_nodal_fields;
real_nodal_fields["damage"] = damage.get();
real_nodal_fields["external_force"] = external_force.get();
real_nodal_fields["internal_force"] = internal_force.get();
return mesh.createNodalField(real_nodal_fields[field_name], group_name);
std::shared_ptr<dumpers::Field> field;
return field;
}
/* -------------------------------------------------------------------------- */
std::shared_ptr<dumpers::Field> PhaseFieldModel::createElementalField(
const std::string & field_name, const std::string & group_name,
bool /*unused*/, Int /*unused*/, ElementKind element_kind) {
if (field_name == "partitions") {
return mesh.createElementalField<Int, dumpers::ElementPartitionField>(
mesh.getConnectivities(), group_name, this->spatial_dimension,
element_kind);
}
std::shared_ptr<dumpers::Field> field;
return field;
}
} // namespace akantu
diff --git a/src/model/phase_field/phase_field_model.hh b/src/model/phase_field/phase_field_model.hh
index 4f1cc8ff3..b1e38b4a2 100644
--- a/src/model/phase_field/phase_field_model.hh
+++ b/src/model/phase_field/phase_field_model.hh
@@ -1,257 +1,261 @@
/**
* Copyright (©) 2018-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 <http://www.gnu.org/licenses/>.
*/
/* -------------------------------------------------------------------------- */
#include "boundary_condition.hh"
#include "constitutive_laws_handler.hh"
#include "fe_engine.hh"
#include "model.hh"
#include "phasefield_selector.hh"
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_PHASE_FIELD_MODEL_HH_
#define AKANTU_PHASE_FIELD_MODEL_HH_
namespace akantu {
class PhaseField;
template <ElementKind kind, class IntegrationOrderFuntor> class IntegratorGauss;
template <ElementKind kind> class ShapeLagrange;
} // namespace akantu
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
class PhaseFieldModel : public ConstitutiveLawsHandler<PhaseField, Model>,
public BoundaryCondition<PhaseFieldModel> {
using Parent = ConstitutiveLawsHandler<PhaseField, Model>;
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
using FEEngineType = FEEngineTemplate<IntegratorGauss, ShapeLagrange>;
PhaseFieldModel(Mesh & mesh, Int dim = _all_dimensions,
const ID & id = "phase_field_model",
std::shared_ptr<DOFManager> dof_manager = nullptr,
ModelType model_type = ModelType::_phase_field_model);
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
protected:
/// generic function to initialize everything ready for explicit dynamics
void initFullImpl(const ModelOptions & options) override;
/// allocate all vectors
void initSolver(TimeStepSolverType /*unused*/,
NonLinearSolverType /*unused*/) override;
/// predictor
void predictor() override;
/// corrector
void corrector() override;
/// compute the heat flux
void assembleResidual() override;
/// get the type of matrix needed
[[nodiscard]] MatrixType getMatrixType(const ID & /*unused*/) const override;
/// callback to assemble a Matrix
void assembleMatrix(const ID & /*unused*/) override;
/// callback to assemble a lumped Matrix
void assembleLumpedMatrix(const ID & /*unused*/) override;
protected:
/* ------------------------------------------------------------------------ */
[[nodiscard]] TimeStepSolverType getDefaultSolverType() const override;
std::tuple<ID, TimeStepSolverType>
getDefaultSolverID(const AnalysisMethod & method) override;
[[nodiscard]] ModelSolverOptions
getDefaultSolverOptions(const TimeStepSolverType & type) const override;
/// set the element_id_by_phasefield and add the elements to the good
/// phasefields
void
assignPhaseFieldToElements(const ElementTypeMapArray<Idx> * filter = nullptr);
/* ------------------------------------------------------------------------ */
/* Methods for static */
/* ------------------------------------------------------------------------ */
public:
/// assembles the phasefield stiffness matrix
virtual void assembleStiffnessMatrix();
/// compute the internal forces
virtual void assembleInternalForces();
// compute the internal forces
void assembleInternalForces(GhostType ghost_type);
// save fields history
void savePreviousState();
/* ------------------------------------------------------------------------ */
/* Methods for dynamic */
/* ------------------------------------------------------------------------ */
public:
/// set the stable timestep
void setTimeStep(Real time_step, const ID & solver_id = "") override;
protected:
/// callback for the solver, this is called at beginning of solve
void beforeSolveStep() override;
/// callback for the solver, this is called at end of solve
void afterSolveStep(bool converged = true) override;
void computeNonLocalContribution(GhostType /*ghost_type*/) override{};
/* ------------------------------------------------------------------------ */
/* Data Accessor inherited members */
/* ------------------------------------------------------------------------ */
public:
[[nodiscard]] Int getNbData(const Array<Element> & elements,
const SynchronizationTag & tag) const override;
void packData(CommunicationBuffer & buffer, const Array<Element> & elements,
const SynchronizationTag & tag) const override;
void unpackData(CommunicationBuffer & buffer, const Array<Element> & elements,
const SynchronizationTag & tag) override;
[[nodiscard]] Int getNbData(const Array<Idx> & indexes,
const SynchronizationTag & tag) const override;
void packData(CommunicationBuffer & buffer, const Array<Idx> & indexes,
const SynchronizationTag & tag) const override;
void unpackData(CommunicationBuffer & buffer, const Array<Idx> & indexes,
const SynchronizationTag & tag) override;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// return the damage array
AKANTU_GET_MACRO_DEREF_PTR(Damage, damage);
AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Damage, damage);
+ /// return the previous damage array
+ AKANTU_GET_MACRO_DEREF_PTR(PreviousDamage, previous_damage);
+ AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(PreviousDamage, previous_damage);
+
/// get the PhaseFieldModel::internal_force vector (internal forces)
AKANTU_GET_MACRO_DEREF_PTR(InternalForce, internal_force);
AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(InternalForce, internal_force);
/// get the PhaseFieldModel::external_force vector (external forces)
AKANTU_GET_MACRO_DEREF_PTR(ExternalForce, external_force);
AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(ExternalForce, external_force);
/// get the PhaseFieldModel::blocked_dofs vector
AKANTU_GET_MACRO_DEREF_PTR(BlockedDOFs, blocked_dofs);
/// get the PhaseFieldModel::blocked_dofs vector
AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(BlockedDOFs, blocked_dofs);
/// Returns the total dissipated energy
[[nodiscard]] virtual Real getEnergy(const ID & energy_id = "dissipated");
/// Compute dissipated energy for an individual element
[[nodiscard]] virtual Real getEnergy(const ID & energy_id,
const Element & element);
/// Compute dissipated energy for an element group
[[nodiscard]] virtual Real getEnergy(const ID & energy_id,
const ID & group_id);
FEEngine & getFEEngineBoundary(const ID & name = "") override;
void
setPhaseFieldSelector(const std::shared_ptr<PhaseFieldSelector> & selector) {
this->setConstitutiveLawSelector(selector);
}
PhaseField & getPhaseField(const ID & id) { return getConstitutiveLaw(id); }
PhaseField & getPhaseField(Idx id) { return getConstitutiveLaw(id); }
[[nodiscard]] const PhaseField & getPhaseField(const ID & id) const {
return getConstitutiveLaw(id);
}
[[nodiscard]] const PhaseField & getPhaseField(Idx id) const {
return getConstitutiveLaw(id);
}
/// give the number of phasefield materials
[[nodiscard]] inline auto getNbPhaseFields() const {
return this->getNbConstitutiveLaws();
}
/* ------------------------------------------------------------------------ */
/* Dumpable Interface */
/* ------------------------------------------------------------------------ */
public:
std::shared_ptr<dumpers::Field>
createNodalFieldReal(const std::string & field_name,
const std::string & group_name,
bool padding_flag) override;
std::shared_ptr<dumpers::Field>
createNodalFieldBool(const std::string & field_name,
const std::string & group_name,
bool padding_flag) override;
std::shared_ptr<dumpers::Field>
createElementalField(const std::string & field_name,
const std::string & group_name, bool padding_flag,
Int spatial_dimension, ElementKind kind) override;
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
/// number of iterations
Int n_iter;
/// damage array
std::unique_ptr<Array<Real>> damage;
/// damage array at the previous time step
std::unique_ptr<Array<Real>> previous_damage;
/// boundary vector
std::unique_ptr<Array<bool>> blocked_dofs;
/// external force vector
std::unique_ptr<Array<Real>> external_force;
/// residuals array
std::unique_ptr<Array<Real>> internal_force;
};
} // namespace akantu
#include "phasefield.hh"
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#include "phase_field_model_inline_impl.hh"
/* -------------------------------------------------------------------------- */
#endif