diff --git a/examples/phase_field/CMakeLists.txt b/examples/phase_field/CMakeLists.txt
index a8e195a16..4e98705d4 100644
--- a/examples/phase_field/CMakeLists.txt
+++ b/examples/phase_field/CMakeLists.txt
@@ -1,32 +1,38 @@
#===============================================================================
# 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 .
#
# @section DESCRIPTION
#
#===============================================================================
add_mesh(square_notch square_notch.geo 2 1)
register_example(phase_field_notch
SOURCES phase_field_notch.cc
DEPENDS square_notch
FILES_TO_COPY material_notch.dat
)
+register_example(phase_field_parallel
+ SOURCES phase_field_parallel.cc
+ DEPENDS square_notch
+ FILES_TO_COPY material_notch.dat
+ )
+
diff --git a/examples/phase_field/material_notch.dat b/examples/phase_field/material_notch.dat
index 5ce8d9c68..a2ace9ea9 100644
--- a/examples/phase_field/material_notch.dat
+++ b/examples/phase_field/material_notch.dat
@@ -1,20 +1,20 @@
model solid_mechanics_model [
material phasefield [
name = plate
rho = 1.
E = 210.0
nu = 0.3
eta = 0
Plane_Stress = false
]
]
model phase_field_model [
phasefield exponential [
name = plate
- l0 = 0.0075
+ l0 = 0.2
gc = 2.7e-3
E = 210.0
nu = 0.3
]
-]
\ No newline at end of file
+]
diff --git a/examples/phase_field/phase_field_notch.cc b/examples/phase_field/phase_field_notch.cc
index de2019ab3..9a9a6a88d 100644
--- a/examples/phase_field/phase_field_notch.cc
+++ b/examples/phase_field/phase_field_notch.cc
@@ -1,109 +1,128 @@
/**
* 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 .
*/
/* -------------------------------------------------------------------------- */
#include "coupler_solid_phasefield.hh"
+#include "group_manager.hh"
#include "non_linear_solver.hh"
+#include "phase_field_element_filter.hh"
#include "phase_field_model.hh"
#include "solid_mechanics_model.hh"
/* -------------------------------------------------------------------------- */
#include
#include
#include
/* -------------------------------------------------------------------------- */
using namespace akantu;
using clk = std::chrono::high_resolution_clock;
using second = std::chrono::duration;
using millisecond = std::chrono::duration;
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");
CouplerSolidPhaseField coupler(mesh);
auto & model = coupler.getSolidMechanicsModel();
auto & phase = coupler.getPhaseFieldModel();
model.initFull(_analysis_method = _static);
auto && mat_selector =
std::make_shared>("physical_names",
model);
model.setMaterialSelector(mat_selector);
auto && selector = std::make_shared>(
"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");
model.addDumpField("stress");
model.addDumpField("grad_u");
model.addDumpFieldVector("displacement");
model.addDumpField("damage");
model.dump();
- auto nbSteps = 1500;
- Real increment = 1e-5;
+ Int nbSteps = 1000;
+ Real increment = 6e-6;
+ Int nb_staggered_steps = 5;
auto start_time = clk::now();
for (Int s = 1; s < nbSteps; ++s) {
if (s >= 500) {
- increment = 1.e-6;
+ increment = 2e-6;
+ nb_staggered_steps = 10;
}
- if (s % 10 == 0) {
+ if (s % 200 == 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");
- coupler.solve();
+ for (Idx i = 0; i < nb_staggered_steps; ++i) {
+ coupler.solve();
+ }
+
+ auto energy = phase.getEnergy();
if (s % 100 == 0) {
model.dump();
}
}
+ Real damage_limit = 0.15;
+ auto global_nb_clusters = mesh.createClusters(
+ spatial_dimension, "crack", PhaseFieldElementFilter(phase, damage_limit));
+
+ auto nb_fragment = mesh.getNbElementGroups(spatial_dimension);
+
+ model.dumpGroup("crack_0");
+
+ std::cout << std::endl;
+ std::cout << "Nb clusters: " << global_nb_clusters << std::endl;
+ std::cout << "Nb fragments: " << nb_fragment << std::endl;
+
finalize();
return EXIT_SUCCESS;
}
diff --git a/examples/phase_field/phase_field_notch.cc b/examples/phase_field/phase_field_parallel.cc
similarity index 75%
copy from examples/phase_field/phase_field_notch.cc
copy to examples/phase_field/phase_field_parallel.cc
index de2019ab3..0d8278f55 100644
--- a/examples/phase_field/phase_field_notch.cc
+++ b/examples/phase_field/phase_field_parallel.cc
@@ -1,109 +1,135 @@
/**
- * Copyright (©) 2018-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
- * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
+ * @file phase_field_parallel.cc
+ *
+ * @author Mohit Pundir
+ *
+ * @date creation: Mon May 09 2022
+ * @date last modification: Mon May 09 2022
+ *
+ * @brief Example of phase field model in parallel
*
- * This file is part of Akantu
+ *
+ * @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 .
+ *
*/
/* -------------------------------------------------------------------------- */
+#include "communicator.hh"
#include "coupler_solid_phasefield.hh"
+#include "group_manager.hh"
#include "non_linear_solver.hh"
#include "phase_field_model.hh"
#include "solid_mechanics_model.hh"
/* -------------------------------------------------------------------------- */
#include
#include
#include
/* -------------------------------------------------------------------------- */
using namespace akantu;
using clk = std::chrono::high_resolution_clock;
using second = std::chrono::duration;
using millisecond = std::chrono::duration;
-const Int spatial_dimension = 2;
-
-/* -------------------------------------------------------------------------- */
+const UInt 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");
+
+ 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>("physical_names",
model);
model.setMaterialSelector(mat_selector);
auto && selector = std::make_shared>(
"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");
+ model.setBaseName("phase_notch_parallel");
model.addDumpField("stress");
model.addDumpField("grad_u");
model.addDumpFieldVector("displacement");
model.addDumpField("damage");
- model.dump();
+ if (mesh.isDistributed()) {
+ // phase.addDumpField("partitions");
+ }
+ phase.dump();
- auto nbSteps = 1500;
- Real increment = 1e-5;
+ UInt nbSteps = 1000;
+ Real increment = 6e-6;
+ UInt nb_staggered_steps = 5;
auto start_time = clk::now();
- for (Int s = 1; s < nbSteps; ++s) {
+ for (UInt s = 1; s < nbSteps; ++s) {
if (s >= 500) {
- increment = 1.e-6;
+ increment = 2e-6;
+ nb_staggered_steps = 10;
}
- if (s % 10 == 0) {
+ 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");
- coupler.solve();
+ for (UInt i = 0; i < nb_staggered_steps; ++i) {
+ coupler.solve();
+ }
if (s % 100 == 0) {
model.dump();
}
}
finalize();
return EXIT_SUCCESS;
}
diff --git a/examples/phase_field/square_notch.geo b/examples/phase_field/square_notch.geo
index d9ab8b202..b3ac5aafa 100644
--- a/examples/phase_field/square_notch.geo
+++ b/examples/phase_field/square_notch.geo
@@ -1,30 +1,30 @@
-element_size = 0.050;
+element_size = 0.1;
Point(1) = {0.5, 0.5, 0, element_size};
Point(2) = {-0.5, 0.5, 0, element_size};
Point(3) = {-0.5, -0.5, 0, element_size};
Point(4) = {0.5, -0.5, 0, element_size};
Point(5) = {-0.5, 0.001, 0, element_size};
Point(6) = {0., 0.0, 0, element_size};
Point(7) = {0.5, 0.0, 0, element_size};
Point(8) = {-0.5, -0.001, 0, element_size};
Line(1) = {3, 4};
Line(2) = {4, 7};
Line(3) = {7, 1};
Line(4) = {1, 2};
Line(5) = {2, 5};
Line(6) = {5, 6};
Line(7) = {6, 8};
Line(8) = {8, 3};
Line Loop(1) = {1, 2, 3, 4, 5, 6, 7, 8};
Plane Surface(1) = {1};
Physical Surface("plate") = {1};
Physical Line("bottom") = {1};
Physical Line("right") = {2, 3};
Physical Line("top") = {4};
Physical Line("left") = {5,8};
diff --git a/packages/phase_field.cmake b/packages/phase_field.cmake
index 83c58ae95..a3d60fd97 100644
--- a/packages/phase_field.cmake
+++ b/packages/phase_field.cmake
@@ -1,44 +1,46 @@
#===============================================================================
# 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 .
#
#===============================================================================
package_declare(phase_field DEFAULT ON
DEPENDS model_couplers
DESCRIPTION "Use Phase Field package of Akantu")
package_declare_sources(phase_field
model/phase_field/phasefield.cc
model/phase_field/phasefield.hh
- model/phase_field/phasefield_inline_impl.cc
+ model/phase_field/phasefield_inline_impl.hh
model/phase_field/phasefield_selector.hh
model/phase_field/phasefield_selector_tmpl.hh
model/phase_field/phasefields/phasefield_exponential.hh
model/phase_field/phasefields/phasefield_exponential.cc
model/phase_field/phasefields/phasefield_exponential_inline_impl.hh
model/phase_field/phase_field_model.cc
model/phase_field/phase_field_model.hh
- model/phase_field/phase_field_model_inline_impl.cc
+ model/phase_field/phase_field_model_inline_impl.hh
model/model_couplers/coupler_solid_phasefield.hh
model/model_couplers/coupler_solid_phasefield.cc
+
+ model/phase_field/phase_field_element_filter.hh
)
diff --git a/python/py_phase_field_model.cc b/python/py_phase_field_model.cc
index 79ce59e1d..6b7071fb9 100644
--- a/python/py_phase_field_model.cc
+++ b/python/py_phase_field_model.cc
@@ -1,132 +1,141 @@
/**
* Copyright (©) 2019-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 "py_aka_array.hh"
/* -------------------------------------------------------------------------- */
#include
#include
#include
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
namespace py = pybind11;
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
#define def_function_nocopy(func_name) \
def( \
#func_name, \
[](PhaseFieldModel & self) -> decltype(auto) { \
return self.func_name(); \
}, \
py::return_value_policy::reference)
#define def_function(func_name) \
def(#func_name, [](PhaseFieldModel & self) -> decltype(auto) { \
return self.func_name(); \
})
/* -------------------------------------------------------------------------- */
void register_phase_field_model(py::module & mod) {
py::class_(mod, "PhaseFieldModelOptions")
.def(py::init(), py::arg("analysis_method") = _static);
py::class_(mod, "PhaseFieldModel",
py::multiple_inheritance())
- .def(py::init(),
+ .def(py::init,
+ const ModelType>(),
py::arg("mesh"), py::arg("spatial_dimension") = _all_dimensions,
py::arg("id") = "phase_field_model",
+ py::arg("dof_manager") = nullptr,
py::arg("model_type") = ModelType::_phase_field_model)
.def(
"initFull",
[](PhaseFieldModel & self, const PhaseFieldModelOptions & options) {
self.initFull(options);
},
py::arg("_analysis_method") = PhaseFieldModelOptions())
.def(
"initFull",
[](PhaseFieldModel & self, const AnalysisMethod & analysis_method) {
self.initFull(_analysis_method = analysis_method);
},
py::arg("_analysis_method"))
.def("applyDirichletBC",
[](PhaseFieldModel & /*self*/) {
PyErr_WarnEx(
PyExc_DeprecationWarning,
"applyDirichletBC() is deprecated, use applyBC instead", 1);
})
.def("applyBC",
[](PhaseFieldModel & self, BC::Dirichlet::DirichletFunctor & func,
const std::string & element_group) {
self.applyBC(func, element_group);
})
.def("applyBC",
[](PhaseFieldModel & self, BC::Neumann::NeumannFunctor & func,
const std::string & element_group) {
self.applyBC(func, element_group);
})
.def("setTimeStep", &PhaseFieldModel::setTimeStep, py::arg("time_step"),
py::arg("solver_id") = "")
+ .def("getEnergy", [](PhaseFieldModel & self) { return self.getEnergy(); })
+ .def(
+ "getEnergy",
+ [](PhaseFieldModel & self, const std::string & group_id) {
+ return self.getEnergy(group_id);
+ },
+ py::arg("group_id"))
.def_function(assembleStiffnessMatrix)
.def_function(assembleInternalForces)
.def_function_nocopy(getDamage)
.def_function_nocopy(getInternalForce)
.def_function_nocopy(getBlockedDOFs)
.def_function_nocopy(getMesh)
.def(
"getPhaseField",
[](PhaseFieldModel & self, Idx phase_field_id) -> decltype(auto) {
return self.getPhaseField(phase_field_id);
},
py::arg("phase_field_id"), py::return_value_policy::reference)
.def(
"getPhaseField",
[](PhaseFieldModel & self,
const ID & phase_field_name) -> decltype(auto) {
return self.getPhaseField(phase_field_name);
},
py::arg("phase_field_name"), py::return_value_policy::reference)
.def("getPhaseFieldIndex", &PhaseFieldModel::getPhaseFieldIndex)
.def("setPhaseFieldSelector", &PhaseFieldModel::setPhaseFieldSelector);
}
void register_phase_field_coupler(py::module & mod) {
py::class_(mod, "CouplerSolidPhaseField")
.def(py::init(),
py::arg("mesh"), py::arg("spatial_dimension") = _all_dimensions,
py::arg("id") = "coupler_solid_phasefield",
py::arg("model_type") = ModelType::_coupler_solid_phasefield)
.def("solve",
[](CouplerSolidPhaseField & self, const ID & solid_solver_id,
const ID & phase_solver_id) {
self.solve(solid_solver_id, phase_solver_id);
})
.def("getSolidMechanicsModel",
&CouplerSolidPhaseField::getSolidMechanicsModel,
py::return_value_policy::reference)
.def("getPhaseFieldModel", &CouplerSolidPhaseField::getPhaseFieldModel,
py::return_value_policy::reference);
}
} // namespace akantu
diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh
index 933696fd1..e7bc42ac7 100644
--- a/src/common/aka_common.hh
+++ b/src/common/aka_common.hh
@@ -1,704 +1,700 @@
/**
* 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 .
*/
/* -------------------------------------------------------------------------- */
#ifndef AKANTU_COMMON_HH_
#define AKANTU_COMMON_HH_
#include "aka_compatibilty_with_cpp_standard.hh"
/* -------------------------------------------------------------------------- */
#if defined(WIN32)
#define __attribute__(x)
#endif
/* -------------------------------------------------------------------------- */
#include "aka_config.hh"
#include "aka_error.hh"
#include "aka_safe_enum.hh"
/* -------------------------------------------------------------------------- */
#include
#include
#include
#include
#include
#include
#include
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
/* Constants */
/* -------------------------------------------------------------------------- */
namespace {
[[gnu::unused]] constexpr Int _all_dimensions{
std::numeric_limits::max()};
#ifdef AKANTU_NDEBUG
[[gnu::unused]] constexpr Real REAL_INIT_VALUE{0.};
#else
[[gnu::unused]] constexpr Real REAL_INIT_VALUE{
std::numeric_limits::quiet_NaN()};
#endif
} // namespace
using dim_1_t = std::integral_constant;
using dim_2_t = std::integral_constant;
using dim_3_t = std::integral_constant;
using AllSpatialDimensions = std::tuple;
/* --------------------------------------------------------------------------
*/
/* Common types */
/* --------------------------------------------------------------------------
*/
using ID = std::string;
} // namespace akantu
/* -------------------------------------------------------------------------- */
#include "aka_enum_macros.hh"
/* -------------------------------------------------------------------------- */
#include "aka_element_classes_info.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
/* Mesh/FEM/Model types */
/* -------------------------------------------------------------------------- */
/// small help to use names for directions
enum SpatialDirection { _x = 0, _y = 1, _z = 2 };
/// enum MeshIOType type of mesh reader/writer
enum MeshIOType {
_miot_auto, ///< Auto guess of the reader to use based on the extension
_miot_gmsh, ///< Gmsh files
_miot_gmsh_struct, ///< Gsmh reader with reintpretation of elements has
/// structures elements
_miot_diana, ///< TNO Diana mesh format
_miot_abaqus ///< Abaqus mesh format
};
/// enum MeshEventHandlerPriority defines relative order of execution of
/// events
enum EventHandlerPriority {
_ehp_highest = 0,
_ehp_mesh = 5,
_ehp_fe_engine = 9,
_ehp_synchronizer = 10,
_ehp_dof_manager = 20,
_ehp_model = 94,
_ehp_non_local_manager = 100,
_ehp_lowest = 100
};
#if !defined(DOXYGEN)
// clang-format off
#define AKANTU_MODEL_TYPES \
(model) \
(solid_mechanics_model) \
(solid_mechanics_model_cohesive) \
(heat_transfer_model) \
(structural_mechanics_model) \
(embedded_model) \
(contact_mechanics_model) \
(coupler_solid_contact) \
(coupler_solid_cohesive_contact) \
(phase_field_model) \
(coupler_solid_phasefield)
// clang-format on
/// enum ModelType defines which type of physics is solved
AKANTU_CLASS_ENUM_DECLARE(ModelType, AKANTU_MODEL_TYPES)
AKANTU_CLASS_ENUM_OUTPUT_STREAM(ModelType, AKANTU_MODEL_TYPES)
AKANTU_CLASS_ENUM_INPUT_STREAM(ModelType, AKANTU_MODEL_TYPES)
#else
enum class ModelType {
model,
solid_mechanics_model,
solid_mechanics_model_cohesive,
heat_transfer_model,
structural_mechanics_model,
embedded_model,
};
#endif
/// enum AnalysisMethod type of solving method used to solve the equation of
/// motion
enum AnalysisMethod {
_static = 0,
_implicit_dynamic = 1,
_explicit_lumped_mass = 2,
_explicit_lumped_capacity = 2,
_explicit_consistent_mass = 3,
_explicit_contact = 4,
_implicit_contact = 5
};
/// enum DOFSupportType defines which kind of dof that can exists
enum DOFSupportType { _dst_nodal, _dst_generic };
#if !defined(DOXYGEN)
// clang-format off
#define AKANTU_NON_LINEAR_SOLVER_TYPES \
(linear) \
(newton_raphson) \
(newton_raphson_modified) \
(lumped) \
(gmres) \
(bfgs) \
(cg) \
(newton_raphson_contact) \
(auto)
// clang-format on
AKANTU_CLASS_ENUM_DECLARE(NonLinearSolverType, AKANTU_NON_LINEAR_SOLVER_TYPES)
AKANTU_CLASS_ENUM_OUTPUT_STREAM(NonLinearSolverType,
AKANTU_NON_LINEAR_SOLVER_TYPES)
AKANTU_CLASS_ENUM_INPUT_STREAM(NonLinearSolverType,
AKANTU_NON_LINEAR_SOLVER_TYPES)
#else
/// Type of non linear resolution available in akantu
enum class NonLinearSolverType {
_linear, ///< No non linear convergence loop
_newton_raphson, ///< Regular Newton-Raphson
_newton_raphson_modified, ///< Newton-Raphson with initial tangent
_lumped, ///< Case of lumped mass or equivalent matrix
_gmres,
_bfgs,
_cg,
_newton_raphson_contact, ///< Regular Newton-Raphson modified
/// for contact problem
_auto, ///< This will take a default value that make sense in case of
/// model::getNewSolver
};
#endif
#if !defined(DOXYGEN)
// clang-format off
#define AKANTU_TIME_STEP_SOLVER_TYPE \
(static) \
(dynamic) \
(dynamic_lumped) \
(not_defined)
// clang-format on
AKANTU_CLASS_ENUM_DECLARE(TimeStepSolverType, AKANTU_TIME_STEP_SOLVER_TYPE)
AKANTU_CLASS_ENUM_OUTPUT_STREAM(TimeStepSolverType,
AKANTU_TIME_STEP_SOLVER_TYPE)
AKANTU_CLASS_ENUM_INPUT_STREAM(TimeStepSolverType, AKANTU_TIME_STEP_SOLVER_TYPE)
#else
/// Type of time stepping solver
enum class TimeStepSolverType {
_static, ///< Static solution
_dynamic, ///< Dynamic solver
_dynamic_lumped, ///< Dynamic solver with lumped mass
_not_defined, ///< For not defined cases
};
#endif
#if !defined(DOXYGEN)
// clang-format off
#define AKANTU_INTEGRATION_SCHEME_TYPE \
(pseudo_time) \
(forward_euler) \
(trapezoidal_rule_1) \
(backward_euler) \
(central_difference) \
(fox_goodwin) \
(trapezoidal_rule_2) \
(linear_acceleration) \
(newmark_beta) \
(generalized_trapezoidal)
// clang-format on
AKANTU_CLASS_ENUM_DECLARE(IntegrationSchemeType, AKANTU_INTEGRATION_SCHEME_TYPE)
AKANTU_CLASS_ENUM_OUTPUT_STREAM(IntegrationSchemeType,
AKANTU_INTEGRATION_SCHEME_TYPE)
AKANTU_CLASS_ENUM_INPUT_STREAM(IntegrationSchemeType,
AKANTU_INTEGRATION_SCHEME_TYPE)
#else
/// Type of integration scheme
enum class IntegrationSchemeType {
_pseudo_time, ///< Pseudo Time
_forward_euler, ///< GeneralizedTrapezoidal(0)
_trapezoidal_rule_1, ///< GeneralizedTrapezoidal(1/2)
_backward_euler, ///< GeneralizedTrapezoidal(1)
_central_difference, ///< NewmarkBeta(0, 1/2)
_fox_goodwin, ///< NewmarkBeta(1/6, 1/2)
_trapezoidal_rule_2, ///< NewmarkBeta(1/2, 1/2)
_linear_acceleration, ///< NewmarkBeta(1/3, 1/2)
_newmark_beta, ///< generic NewmarkBeta with user defined
/// alpha and beta
_generalized_trapezoidal ///< generic GeneralizedTrapezoidal with user
/// defined alpha
};
#endif
#if !defined(DOXYGEN)
// clang-format off
#define AKANTU_SOLVE_CONVERGENCE_CRITERIA \
(residual) \
(solution) \
(residual_mass_wgh)
// clang-format on
AKANTU_CLASS_ENUM_DECLARE(SolveConvergenceCriteria,
AKANTU_SOLVE_CONVERGENCE_CRITERIA)
AKANTU_CLASS_ENUM_OUTPUT_STREAM(SolveConvergenceCriteria,
AKANTU_SOLVE_CONVERGENCE_CRITERIA)
AKANTU_CLASS_ENUM_INPUT_STREAM(SolveConvergenceCriteria,
AKANTU_SOLVE_CONVERGENCE_CRITERIA)
#else
/// enum SolveConvergenceCriteria different convergence criteria
enum class SolveConvergenceCriteria {
_residual, ///< Use residual to test the convergence
_solution, ///< Use solution to test the convergence
_residual_mass_wgh ///< Use residual weighted by inv. nodal mass to
///< testb
};
#endif
/// enum CohesiveMethod type of insertion of cohesive elements
enum CohesiveMethod { _intrinsic, _extrinsic };
/// @enum MatrixType type of sparse matrix used
enum MatrixType { _unsymmetric, _symmetric, _mt_not_defined };
/// @enum Type of contact detection
enum DetectionType { _explicit, _implicit };
#if !defined(DOXYGEN)
// clang-format off
#define AKANTU_CONTACT_STATE \
(no_contact) \
(stick) \
(slip)
// clang-format on
AKANTU_CLASS_ENUM_DECLARE(ContactState, AKANTU_CONTACT_STATE)
AKANTU_CLASS_ENUM_OUTPUT_STREAM(ContactState, AKANTU_CONTACT_STATE)
AKANTU_CLASS_ENUM_INPUT_STREAM(ContactState, AKANTU_CONTACT_STATE)
#else
/// @enum no contact or stick or slip state
enum class ContactState {
_no_contact = 0,
_stick = 1,
_slip = 2,
};
#endif
/* -------------------------------------------------------------------------- */
/* Ghosts handling */
/* -------------------------------------------------------------------------- */
/// @enum CommunicatorType type of communication method to use
enum CommunicatorType { _communicator_mpi, _communicator_dummy };
#if !defined(DOXYGEN)
// clang-format off
#define AKANTU_SYNCHRONIZATION_TAG \
(whatever) \
(update) \
(ask_nodes) \
(size) \
(smm_mass) \
(smm_for_gradu) \
(smm_boundary) \
(smm_uv) \
(smm_res) \
(smm_init_mat) \
(smm_stress) \
+ (smm_gradu) \
(smmc_facets) \
(smmc_facets_conn) \
(smmc_facets_stress) \
(smmc_damage) \
(giu_global_conn) \
(ce_groups) \
(ce_insertion_order) \
(gm_clusters) \
(htm_temperature) \
(htm_gradient_temperature) \
(htm_phi) \
(htm_gradient_phi) \
(pfm_damage) \
- (pfm_driving) \
- (pfm_history) \
- (pfm_energy) \
(csp_damage) \
(csp_strain) \
(mnl_for_average) \
(mnl_weight) \
(nh_criterion) \
(test) \
(user_1) \
(user_2) \
(material_id) \
+ (phasefield_id) \
(for_dump) \
(cf_nodal) \
(cf_incr) \
(solver_solution)
// clang-format on
AKANTU_CLASS_ENUM_DECLARE(SynchronizationTag, AKANTU_SYNCHRONIZATION_TAG)
AKANTU_CLASS_ENUM_OUTPUT_STREAM(SynchronizationTag, AKANTU_SYNCHRONIZATION_TAG)
#else
/// @enum SynchronizationTag type of synchronizations
enum class SynchronizationTag {
//--- Generic tags ---
_whatever,
_update,
_ask_nodes,
_size,
//--- SolidMechanicsModel tags ---
_smm_mass, ///< synchronization of the SolidMechanicsModel.mass
_smm_for_gradu, ///< synchronization of the
/// SolidMechanicsModel.displacement
_smm_boundary, ///< synchronization of the boundary, forces, velocities
/// and displacement
_smm_uv, ///< synchronization of the nodal velocities and displacement
_smm_res, ///< synchronization of the nodal residual
_smm_init_mat, ///< synchronization of the data to initialize materials
_smm_stress, ///< synchronization of the stresses to compute the
///< internal
/// forces
+ _smm_gradu, ///< synchronization of the gradu to compute the
+ ///< strain
_smmc_facets, ///< synchronization of facet data to setup facet synch
_smmc_facets_conn, ///< synchronization of facet global connectivity
_smmc_facets_stress, ///< synchronization of facets' stress to setup
///< facet
/// synch
_smmc_damage, ///< synchronization of damage
// --- GlobalIdsUpdater tags ---
_giu_global_conn, ///< synchronization of global connectivities
// --- CohesiveElementInserter tags ---
_ce_groups, ///< synchronization of cohesive element insertion depending
/// on facet groups
_ce_insertion_order, ///< synchronization of the order of insertion of
/// cohesive elements
// --- GroupManager tags ---
_gm_clusters, ///< synchronization of clusters
// --- HeatTransfer tags ---
_htm_temperature, ///< synchronization of the nodal temperature
_htm_gradient_temperature, ///< synchronization of the element gradient
/// temperature
// --- PhaseFieldModel tags ---
_pfm_damage, ///< synchronization of the nodal damage
- _pfm_driving, ///< synchronization of the driving forces to
- /// compute the internal
- _pfm_history, ///< synchronization of the damage history to
- /// compute the internal
- _pfm_energy, ///< synchronization of the damage energy
- /// density to compute the internal
-
+
// --- CouplerSolidPhaseField tags ---
_csp_damage, ///< synchronization of the damage from phase
/// model to solid model
_csp_strain, ///< synchronization of the strain from solid
/// model to phase model
// --- LevelSet tags ---
_htm_phi, ///< synchronization of the nodal level set value phi
_htm_gradient_phi, ///< synchronization of the element gradient phi
//--- Material non local ---
_mnl_for_average, ///< synchronization of data to average in non local
/// material
_mnl_weight, ///< synchronization of data for the weight computations
// --- NeighborhoodSynchronization tags ---
_nh_criterion,
// --- General tags ---
_test, ///< Test tag
_user_1, ///< tag for user simulations
_user_2, ///< tag for user simulations
_material_id, ///< synchronization of the material ids
+ _phasefield_id, ///< synchronization of the phasefield ids
_for_dump, ///< everything that needs to be synch before dump
// --- Contact & Friction ---
_cf_nodal, ///< synchronization of disp, velo, and current position
_cf_incr, ///< synchronization of increment
// --- Solver tags ---
_solver_solution ///< synchronization of the solution obained with the
/// PETSc solver
};
#endif
/// @enum GhostType type of ghost
enum GhostType {
_not_ghost = 0,
_ghost = 1,
_casper // not used but a real cute ghost
};
/// Define the flag that can be set to a node
enum class NodeFlag : std::uint8_t {
_normal = 0x00,
_distributed = 0x01,
_master = 0x03,
_slave = 0x05,
_pure_ghost = 0x09,
_shared_mask = 0x0F,
_periodic = 0x10,
_periodic_master = 0x30,
_periodic_slave = 0x50,
_periodic_mask = 0xF0,
_local_master_mask = 0xCC, // ~(_master & _periodic_mask)
};
inline NodeFlag operator&(const NodeFlag & a, const NodeFlag & b) {
using under = std::underlying_type_t;
return NodeFlag(under(a) & under(b));
}
inline NodeFlag operator|(const NodeFlag & a, const NodeFlag & b) {
using under = std::underlying_type_t;
return NodeFlag(under(a) | under(b));
}
inline NodeFlag & operator|=(NodeFlag & a, const NodeFlag & b) {
a = a | b;
return a;
}
inline NodeFlag & operator&=(NodeFlag & a, const NodeFlag & b) {
a = a & b;
return a;
}
inline NodeFlag operator~(const NodeFlag & a) {
using under = std::underlying_type_t;
return NodeFlag(~under(a));
}
std::ostream & operator<<(std::ostream & stream, NodeFlag flag);
} // namespace akantu
AKANTU_ENUM_HASH(GhostType)
namespace akantu {
/* -------------------------------------------------------------------------- */
struct GhostType_def {
using type = GhostType;
static const type _begin_ = _not_ghost;
static const type _end_ = _casper;
};
using ghost_type_t = safe_enum;
namespace {
constexpr ghost_type_t ghost_types{_casper};
}
/// standard output stream operator for GhostType
// inline std::ostream & operator<<(std::ostream & stream, GhostType type);
/* -------------------------------------------------------------------------- */
/* Global defines */
/* -------------------------------------------------------------------------- */
#define AKANTU_MIN_ALLOCATION 2000
#define AKANTU_INDENT ' '
#define AKANTU_INCLUDE_INLINE_IMPL
/* -------------------------------------------------------------------------- */
#define AKANTU_SET_MACRO(name, variable, type) \
inline void set##name(type variable) { this->variable = variable; }
#define AKANTU_GET_MACRO(name, variable, type) \
inline auto get##name() const->type { return variable; }
#define AKANTU_GET_MACRO_AUTO(name, variable) \
inline decltype(auto) get##name() const { return (variable); }
#define AKANTU_GET_MACRO_AUTO_NOT_CONST(name, variable) \
inline decltype(auto) get##name() { return (variable); }
#define AKANTU_GET_MACRO_NOT_CONST(name, variable, type) \
inline auto get##name()->type { return variable; }
#define AKANTU_GET_MACRO_DEREF_PTR(name, ptr) \
inline const auto & get##name() const { \
if (not(ptr)) { \
AKANTU_EXCEPTION("The member " << #ptr << " is not initialized"); \
} \
return (*(ptr)); \
}
#define AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(name, ptr) \
inline decltype(auto) get##name() { \
if (not(ptr)) { \
AKANTU_EXCEPTION("The member " << #ptr << " is not initialized"); \
} \
return (*(ptr)); \
}
#define AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, support, con) \
inline auto get##name(const support & el_type, \
GhostType ghost_type = _not_ghost) \
con->con Array & { \
return variable(el_type, ghost_type); \
} // NOLINT
#define AKANTU_GET_MACRO_BY_ELEMENT_TYPE(name, variable, type) \
AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, )
#define AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(name, variable, type) \
AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, const)
#define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE(name, variable, type) \
AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, )
#define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE_CONST(name, variable, type) \
AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, const)
/* -------------------------------------------------------------------------- */
/// initialize the static part of akantu
void initialize(int & argc, char **& argv);
/// initialize the static part of akantu and read the global input_file
void initialize(const std::string & input_file, int & argc, char **& argv);
/* -------------------------------------------------------------------------- */
/// finilize correctly akantu and clean the memory
void finalize();
/* -------------------------------------------------------------------------- */
/// Read an new input file
void readInputFile(const std::string & input_file);
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/* string manipulation */
/* -------------------------------------------------------------------------- */
inline auto to_lower(const std::string & str) -> std::string;
/* -------------------------------------------------------------------------- */
inline auto trim(const std::string & to_trim) -> std::string;
inline auto trim(const std::string & to_trim, char c) -> std::string;
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/// give a string representation of the a human readable size in bit
template auto printMemorySize(UInt size) -> std::string;
/* -------------------------------------------------------------------------- */
} // namespace akantu
/* -------------------------------------------------------------------------- */
/* Type traits */
/* -------------------------------------------------------------------------- */
namespace aka {
/* ------------------------------------------------------------------------ */
template using is_scalar = std::is_arithmetic;
/* ------------------------------------------------------------------------ */
template > * = nullptr>
auto is_of_type(T && t) -> bool {
return (dynamic_cast>,
std::add_const_t, R>>>(&t) != nullptr);
}
/* -------------------------------------------------------------------------- */
template
auto is_of_type(std::unique_ptr & t) -> bool {
return (dynamic_cast, std::add_const_t, R>>>(
t.get()) != nullptr);
}
/* -------------------------------------------------------------------------- */
template
decltype(auto) as_type(const std::shared_ptr & t) {
return std::dynamic_pointer_cast(t);
}
/* ------------------------------------------------------------------------ */
template > * = nullptr>
decltype(auto) as_type(T && t) {
static_assert(
disjunction<
std::is_base_of, std::decay_t>, // down-cast
std::is_base_of, std::decay_t> // up-cast
>::value,
"Type T and R are not valid for a as_type conversion");
return dynamic_cast>::value,
std::add_const_t, R>>>(t);
}
/* -------------------------------------------------------------------------- */
template ::value> * = nullptr>
decltype(auto) as_type(T && t) {
return &as_type(*t);
}
template inline constexpr auto decay_v = std::decay_t::value;
} // namespace aka
#include "aka_common_inline_impl.hh"
#include "aka_fwd.hh"
namespace akantu {
/// get access to the internal argument parser
cppargparse::ArgumentParser & getStaticArgumentParser();
/// get access to the internal input file parser
Parser & getStaticParser();
/// get access to the user part of the internal input file parser
const ParserSection & getUserParser();
#define AKANTU_CURRENT_FUNCTION \
(std::string(__func__) + "():" + std::to_string(__LINE__))
} // namespace akantu
/* -------------------------------------------------------------------------- */
#if AKANTU_INTEGER_SIZE == 4
#define AKANTU_HASH_COMBINE_MAGIC_NUMBER 0x9e3779b9
#elif AKANTU_INTEGER_SIZE == 8
#define AKANTU_HASH_COMBINE_MAGIC_NUMBER 0x9e3779b97f4a7c13LL
#endif
namespace std {
/**
* Hashing function for pairs based on hash_combine from boost The magic
* number is coming from the golden number @f[\phi = \frac{1 + \sqrt5}{2}@f]
* @f[\frac{2^32}{\phi} = 0x9e3779b9@f]
* http://stackoverflow.com/questions/4948780/magic-number-in-boosthash-combine
* http://burtleburtle.net/bob/hash/doobs.html
*/
template struct hash> {
hash() = default;
auto operator()(const std::pair & p) const -> std::size_t {
size_t seed = ah(p.first);
// NOLINTNEXTLINE(cppcoreguidelines-avoid-magic-numbers)
return bh(p.second) + AKANTU_HASH_COMBINE_MAGIC_NUMBER + (seed << 6) +
(seed >> 2);
}
private:
const hash ah{};
const hash bh{};
};
} // namespace std
#endif // AKANTU_COMMON_HH_
diff --git a/src/model/model_couplers/coupler_solid_phasefield.cc b/src/model/model_couplers/coupler_solid_phasefield.cc
index 87e23f8f4..379e24612 100644
--- a/src/model/model_couplers/coupler_solid_phasefield.cc
+++ b/src/model/model_couplers/coupler_solid_phasefield.cc
@@ -1,567 +1,565 @@
/**
* Copyright (©) 2019-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 "coupler_solid_phasefield.hh"
#include "dumpable_inline_impl.hh"
#include "element_synchronizer.hh"
#include "integrator_gauss.hh"
#include "shape_lagrange.hh"
/* -------------------------------------------------------------------------- */
#include "dumper_iohelper_paraview.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
CouplerSolidPhaseField::CouplerSolidPhaseField(Mesh & mesh, Int dim,
const ID & id,
const ModelType model_type)
: Model(mesh, model_type, dim, id) {
this->registerFEEngineObject("CouplerSolidPhaseField", mesh,
Model::spatial_dimension);
this->mesh.registerDumper("coupler_solid_phasefield", id,
true);
this->mesh.addDumpMeshToDumper("coupler_solid_phasefield", mesh,
Model::spatial_dimension, _not_ghost,
_ek_regular);
this->registerDataAccessor(*this);
solid = new SolidMechanicsModel(mesh, Model::spatial_dimension,
"solid_mechanics_model");
phase =
new PhaseFieldModel(mesh, Model::spatial_dimension, "phase_field_model");
if (this->mesh.isDistributed()) {
auto & synchronizer = this->mesh.getElementSynchronizer();
this->registerSynchronizer(synchronizer, SynchronizationTag::_csp_damage);
this->registerSynchronizer(synchronizer, SynchronizationTag::_csp_strain);
}
}
/* -------------------------------------------------------------------------- */
CouplerSolidPhaseField::~CouplerSolidPhaseField() = default;
/* -------------------------------------------------------------------------- */
void CouplerSolidPhaseField::initFullImpl(const ModelOptions & options) {
Model::initFullImpl(options);
this->initBC(*this, *displacement, *displacement_increment, *external_force);
solid->initFull(_analysis_method = this->method);
phase->initFull(_analysis_method = this->method);
}
/* -------------------------------------------------------------------------- */
void CouplerSolidPhaseField::initModel() {
getFEEngine().initShapeFunctions(_not_ghost);
getFEEngine().initShapeFunctions(_ghost);
}
/* -------------------------------------------------------------------------- */
FEEngine & CouplerSolidPhaseField::getFEEngineBoundary(const ID & name) {
return dynamic_cast(
getFEEngineClassBoundary(name));
}
/* -------------------------------------------------------------------------- */
void CouplerSolidPhaseField::initSolver(
TimeStepSolverType time_step_solver_type,
NonLinearSolverType non_linear_solver_type) {
auto & solid_model_solver = aka::as_type(*solid);
solid_model_solver.initSolver(time_step_solver_type, non_linear_solver_type);
auto & phase_model_solver = aka::as_type(*phase);
phase_model_solver.initSolver(time_step_solver_type, non_linear_solver_type);
}
/* -------------------------------------------------------------------------- */
std::tuple
CouplerSolidPhaseField::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);
}
}
/* -------------------------------------------------------------------------- */
TimeStepSolverType CouplerSolidPhaseField::getDefaultSolverType() const {
return TimeStepSolverType::_dynamic_lumped;
}
/* -------------------------------------------------------------------------- */
ModelSolverOptions CouplerSolidPhaseField::getDefaultSolverOptions(
const TimeStepSolverType & type) const {
ModelSolverOptions options;
switch (type) {
case TimeStepSolverType::_dynamic_lumped: {
options.non_linear_solver_type = NonLinearSolverType::_lumped;
options.integration_scheme_type["displacement"] =
IntegrationSchemeType::_central_difference;
options.solution_type["displacement"] = IntegrationScheme::_acceleration;
break;
}
case TimeStepSolverType::_dynamic: {
options.non_linear_solver_type = NonLinearSolverType::_lumped;
options.integration_scheme_type["displacement"] =
IntegrationSchemeType::_central_difference;
options.solution_type["displacement"] = IntegrationScheme::_acceleration;
break;
}
case TimeStepSolverType::_static: {
options.non_linear_solver_type = NonLinearSolverType::_linear;
options.integration_scheme_type["displacement"] =
IntegrationSchemeType::_pseudo_time;
options.solution_type["displacement"] = IntegrationScheme::_not_defined;
break;
}
default:
AKANTU_EXCEPTION(type << " is not a valid time step solver type");
break;
}
return options;
}
/* -------------------------------------------------------------------------- */
void CouplerSolidPhaseField::assembleResidual() {
// computes the internal forces
this->assembleInternalForces();
auto & solid_internal_force = solid->getInternalForce();
auto & solid_external_force = solid->getExternalForce();
auto & phasefield_internal_force = phase->getInternalForce();
auto & phasefield_external_force = phase->getExternalForce();
/* ------------------------------------------------------------------------ */
this->getDOFManager().assembleToResidual("displacement", solid_external_force,
1);
this->getDOFManager().assembleToResidual("displacement", solid_internal_force,
1);
this->getDOFManager().assembleToResidual("damage", phasefield_external_force,
1);
this->getDOFManager().assembleToResidual("damage", phasefield_internal_force,
1);
}
/* -------------------------------------------------------------------------- */
void CouplerSolidPhaseField::assembleResidual(const ID & residual_part) {
AKANTU_DEBUG_IN();
auto & solid_internal_force = solid->getInternalForce();
auto & solid_external_force = solid->getExternalForce();
auto & phasefield_internal_force = phase->getInternalForce();
auto & phasefield_external_force = phase->getExternalForce();
if ("external" == residual_part) {
this->getDOFManager().assembleToResidual("displacement",
solid_external_force, 1);
this->getDOFManager().assembleToResidual("displacement",
solid_internal_force, 1);
AKANTU_DEBUG_OUT();
return;
}
if ("internal" == residual_part) {
this->getDOFManager().assembleToResidual("damage",
phasefield_external_force, 1);
this->getDOFManager().assembleToResidual("damage",
phasefield_internal_force, 1);
AKANTU_DEBUG_OUT();
return;
}
AKANTU_CUSTOM_EXCEPTION(
debug::SolverCallbackResidualPartUnknown(residual_part));
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void CouplerSolidPhaseField::predictor() {
auto & solid_model_solver = aka::as_type(*solid);
solid_model_solver.predictor();
auto & phase_model_solver = aka::as_type(*phase);
phase_model_solver.predictor();
}
/* -------------------------------------------------------------------------- */
void CouplerSolidPhaseField::corrector() {
auto & solid_model_solver = aka::as_type(*solid);
solid_model_solver.corrector();
auto & phase_model_solver = aka::as_type(*phase);
phase_model_solver.corrector();
}
/* -------------------------------------------------------------------------- */
MatrixType CouplerSolidPhaseField::getMatrixType(const ID & matrix_id) const {
if (matrix_id == "K") {
return _symmetric;
}
if (matrix_id == "M") {
return _symmetric;
}
return _mt_not_defined;
}
/* -------------------------------------------------------------------------- */
void CouplerSolidPhaseField::assembleMatrix(const ID & matrix_id) {
if (matrix_id == "K") {
this->assembleStiffnessMatrix();
} else if (matrix_id == "M") {
solid->assembleMass();
}
}
/* -------------------------------------------------------------------------- */
void CouplerSolidPhaseField::assembleLumpedMatrix(const ID & matrix_id) {
if (matrix_id == "M") {
solid->assembleMassLumped();
}
}
/* -------------------------------------------------------------------------- */
void CouplerSolidPhaseField::beforeSolveStep() {
auto & solid_solver_callback = aka::as_type(*solid);
solid_solver_callback.beforeSolveStep();
auto & phase_solver_callback = aka::as_type(*phase);
phase_solver_callback.beforeSolveStep();
}
/* -------------------------------------------------------------------------- */
void CouplerSolidPhaseField::afterSolveStep(bool converged) {
auto & solid_solver_callback = aka::as_type(*solid);
solid_solver_callback.afterSolveStep(converged);
auto & phase_solver_callback = aka::as_type(*phase);
phase_solver_callback.afterSolveStep(converged);
}
/* -------------------------------------------------------------------------- */
void CouplerSolidPhaseField::assembleInternalForces() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Assemble the internal forces");
solid->assembleInternalForces();
phase->assembleInternalForces();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void CouplerSolidPhaseField::assembleStiffnessMatrix() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Assemble the new stiffness matrix");
solid->assembleStiffnessMatrix();
phase->assembleStiffnessMatrix();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void CouplerSolidPhaseField::assembleMassLumped() {
solid->assembleMassLumped();
}
/* -------------------------------------------------------------------------- */
void CouplerSolidPhaseField::assembleMass() { solid->assembleMass(); }
/* -------------------------------------------------------------------------- */
void CouplerSolidPhaseField::assembleMassLumped(GhostType ghost_type) {
solid->assembleMassLumped(ghost_type);
}
/* -------------------------------------------------------------------------- */
void CouplerSolidPhaseField::assembleMass(GhostType ghost_type) {
solid->assembleMass(ghost_type);
}
/* ------------------------------------------------------------------------- */
void CouplerSolidPhaseField::computeDamageOnQuadPoints(GhostType ghost_type) {
AKANTU_DEBUG_IN();
auto & fem = phase->getFEEngine();
auto & mesh = phase->getMesh();
auto nb_materials = solid->getNbMaterials();
auto nb_phasefields = phase->getNbPhaseFields();
AKANTU_DEBUG_ASSERT(
nb_phasefields == nb_materials,
"The number of phasefields and materials should be equal");
for (auto index : arange(nb_materials)) {
auto & material = solid->getMaterial(index);
for (auto index2 : arange(nb_phasefields)) {
auto & phasefield = phase->getPhaseField(index2);
if (phasefield.getName() == material.getName()) {
switch (spatial_dimension) {
case 1: {
auto & mat = static_cast &>(material);
auto & damage = mat.getDamage();
for (const auto & type :
mesh.elementTypes(Model::spatial_dimension, ghost_type)) {
auto & damage_on_qpoints_vect = damage(type, ghost_type);
fem.interpolateOnIntegrationPoints(phase->getDamage(),
damage_on_qpoints_vect, 1, type,
ghost_type);
}
break;
}
case 2: {
auto & mat = static_cast &>(material);
auto & damage = mat.getDamage();
for (const auto & type :
mesh.elementTypes(Model::spatial_dimension, ghost_type)) {
auto & damage_on_qpoints_vect = damage(type, ghost_type);
fem.interpolateOnIntegrationPoints(phase->getDamage(),
damage_on_qpoints_vect, 1, type,
ghost_type);
}
break;
}
default:
auto & mat = static_cast &>(material);
auto & damage = mat.getDamage();
for (const auto & type :
mesh.elementTypes(Model::spatial_dimension, ghost_type)) {
auto & damage_on_qpoints_vect = damage(type, ghost_type);
fem.interpolateOnIntegrationPoints(phase->getDamage(),
damage_on_qpoints_vect, 1, type,
ghost_type);
}
break;
}
}
}
}
AKANTU_DEBUG_OUT();
}
/* ------------------------------------------------------------------------- */
void CouplerSolidPhaseField::computeStrainOnQuadPoints(GhostType ghost_type) {
AKANTU_DEBUG_IN();
- auto & mesh = solid->getMesh();
-
- auto nb_materials = solid->getNbMaterials();
- auto nb_phasefields = phase->getNbPhaseFields();
-
- AKANTU_DEBUG_ASSERT(
- nb_phasefields == nb_materials,
- "The number of phasefields and materials should be equal");
-
- for (auto index : arange(nb_materials)) {
- auto & material = solid->getMaterial(index);
-
- for (auto index2 : arange(nb_phasefields)) {
- auto & phasefield = phase->getPhaseField(index2);
+ auto & gradu_internal =
+ solid->flattenInternal("grad_u", _ek_regular, ghost_type);
- if (phasefield.getName() == material.getName()) {
-
- auto & strain_on_qpoints = phasefield.getStrain();
- const auto & gradu_on_qpoints = material.getGradU();
-
- for (const auto & type :
- mesh.elementTypes(spatial_dimension, ghost_type)) {
- auto & strain_on_qpoints_vect = strain_on_qpoints(type, ghost_type);
- const auto & gradu_on_qpoints_vect =
- gradu_on_qpoints(type, ghost_type);
- for (auto && values :
- zip(make_view(strain_on_qpoints_vect, spatial_dimension,
- spatial_dimension),
- make_view(gradu_on_qpoints_vect, spatial_dimension,
- spatial_dimension))) {
- auto & strain = std::get<0>(values);
- const auto & grad_u = std::get<1>(values);
- strain = (grad_u + grad_u.transpose()) / 2.;
- }
- }
-
- break;
- }
+ auto & mesh = solid->getMesh();
+ auto & fem = solid->getFEEngine();
+
+ ElementTypeMapArray strain_tmp("temporary strain on quads");
+ strain_tmp.initialize(
+ fem, _nb_component = spatial_dimension * spatial_dimension,
+ _spatial_dimension = spatial_dimension, _ghost_type = ghost_type);
+
+ for (const auto & type : mesh.elementTypes(spatial_dimension, ghost_type)) {
+ auto & strain_vect = strain_tmp(type, ghost_type);
+ const auto & gradu_vect = gradu_internal(type, ghost_type);
+ for (auto && values :
+ zip(make_view(gradu_vect, spatial_dimension, spatial_dimension),
+ make_view(strain_vect, spatial_dimension, spatial_dimension))) {
+ const auto & grad_u = std::get<0>(values);
+ auto & strain = std::get<1>(values);
+ strain = (grad_u + grad_u.transpose()) / 2.;
}
}
+ phase->inflateInternal("strain", strain_tmp, _ek_regular, ghost_type);
+
AKANTU_DEBUG_OUT();
}
/* ------------------------------------------------------------------------- */
void CouplerSolidPhaseField::solve(const ID & solid_solver_id,
const ID & phase_solver_id) {
solid->solveStep(solid_solver_id);
+
+ solid->synchronize(SynchronizationTag::_smm_gradu);
+
+ AKANTU_DEBUG_INFO("exchange strain for local elements");
this->computeStrainOnQuadPoints(_not_ghost);
+ AKANTU_DEBUG_INFO("exchange strain for ghost elements");
+ this->computeStrainOnQuadPoints(_ghost);
+
phase->solveStep(phase_solver_id);
+
+ phase->synchronize(SynchronizationTag::_pfm_damage);
+
+ AKANTU_DEBUG_INFO("exchange damage for local elements");
this->computeDamageOnQuadPoints(_not_ghost);
+ AKANTU_DEBUG_INFO("exchange damage for ghost elements");
+ this->computeDamageOnQuadPoints(_ghost);
+
solid->assembleInternalForces();
}
/* ------------------------------------------------------------------------- */
bool CouplerSolidPhaseField::checkConvergence(Array & u_new,
Array & u_old,
Array & d_new,
Array & d_old) {
const Array & blocked_dofs = solid->getBlockedDOFs();
Int nb_degree_of_freedom = u_new.size();
auto u_n_it = u_new.begin();
auto u_o_it = u_old.begin();
auto bld_it = blocked_dofs.begin();
Real norm = 0;
for (Int n = 0; n < nb_degree_of_freedom; ++n, ++u_n_it, ++u_o_it, ++bld_it) {
if ((!*bld_it)) {
norm += (*u_n_it - *u_o_it) * (*u_n_it - *u_o_it);
}
}
norm = std::sqrt(norm);
auto d_n_it = d_new.begin();
auto d_o_it = d_old.begin();
nb_degree_of_freedom = d_new.size();
Real norm2 = 0;
for (Int i = 0; i < nb_degree_of_freedom; ++i) {
norm2 += (*d_n_it - *d_o_it);
}
norm2 = std::sqrt(norm2);
Real error = std::max(norm, norm2);
Real tolerance = 1e-8;
return error < tolerance;
}
/* -------------------------------------------------------------------------- */
std::shared_ptr CouplerSolidPhaseField::createElementalField(
const std::string & field_name, const std::string & group_name,
bool padding_flag, Idx spatial_dimension, ElementKind kind) {
return solid->createElementalField(field_name, group_name, padding_flag,
spatial_dimension, kind);
std::shared_ptr field;
return field;
}
/* -------------------------------------------------------------------------- */
std::shared_ptr
CouplerSolidPhaseField::createNodalFieldReal(const std::string & field_name,
const std::string & group_name,
bool padding_flag) {
return solid->createNodalFieldReal(field_name, group_name, padding_flag);
std::shared_ptr field;
return field;
}
/* -------------------------------------------------------------------------- */
std::shared_ptr
CouplerSolidPhaseField::createNodalFieldBool(const std::string & field_name,
const std::string & group_name,
bool padding_flag) {
return solid->createNodalFieldBool(field_name, group_name, padding_flag);
std::shared_ptr field;
return field;
}
/* -----------------------------------------------------------------------*/
void CouplerSolidPhaseField::dump(const std::string & dumper_name) {
solid->onDump();
mesh.dump(dumper_name);
}
/* ------------------------------------------------------------------------*/
void CouplerSolidPhaseField::dump(const std::string & dumper_name, Int step) {
solid->onDump();
mesh.dump(dumper_name, step);
}
/* ----------------------------------------------------------------------- */
void CouplerSolidPhaseField::dump(const std::string & dumper_name, Real time,
Int step) {
solid->onDump();
mesh.dump(dumper_name, time, step);
}
/* -------------------------------------------------------------------------- */
void CouplerSolidPhaseField::dump() {
solid->onDump();
mesh.dump();
}
/* -------------------------------------------------------------------------- */
void CouplerSolidPhaseField::dump(Int step) {
solid->onDump();
mesh.dump(step);
}
/* -------------------------------------------------------------------------- */
void CouplerSolidPhaseField::dump(Real time, Int step) {
solid->onDump();
mesh.dump(time, step);
}
} // namespace akantu
diff --git a/src/model/model_couplers/coupler_solid_phasefield.hh b/src/model/model_couplers/coupler_solid_phasefield.hh
index 44f4c2787..20894e7d9 100644
--- a/src/model/model_couplers/coupler_solid_phasefield.hh
+++ b/src/model/model_couplers/coupler_solid_phasefield.hh
@@ -1,246 +1,247 @@
/**
* Copyright (©) 2019-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 "boundary_condition.hh"
#include "data_accessor.hh"
#include "fe_engine.hh"
#include "material.hh"
#include "material_phasefield.hh"
#include "model.hh"
+#include "phasefield.hh"
#include "phase_field_model.hh"
#include "solid_mechanics_model.hh"
#include "sparse_matrix.hh"
#include "time_step_solver.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_COUPLER_SOLID_PHASEFIELD_HH__
#define __AKANTU_COUPLER_SOLID_PHASEFIELD_HH__
/* ------------------------------------------------------------------------ */
/* Coupling : Solid Mechanics / PhaseField */
/* ------------------------------------------------------------------------ */
namespace akantu {
template
class IntegratorGauss;
template class ShapeLagrange;
class DOFManager;
} // namespace akantu
namespace akantu {
class CouplerSolidPhaseField
: public Model,
public DataAccessor,
public DataAccessor,
public BoundaryCondition {
/* ------------------------------------------------------------------------ */
/* Constructor/Destructors */
/* ------------------------------------------------------------------------ */
using MyFEEngineType = FEEngineTemplate;
public:
CouplerSolidPhaseField(
Mesh & mesh, Int dim = _all_dimensions,
const ID & id = "coupler_solid_phasefield",
ModelType model_type = ModelType::_coupler_solid_phasefield);
~CouplerSolidPhaseField() override;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
protected:
/// initialize the complete model
void initFullImpl(const ModelOptions & options) override;
/// initialize the modelType
void initModel() override;
/// get some default values for derived classes
std::tuple
getDefaultSolverID(const AnalysisMethod & method) override;
/* ------------------------------------------------------------------------ */
/* Solver Interface */
/* ------------------------------------------------------------------------ */
public:
/// assembles the contact stiffness matrix
virtual void assembleStiffnessMatrix();
/// assembles the contant internal forces
virtual void assembleInternalForces();
public:
/// computes damage on quad points for solid mechanics model from
/// damage array from phasefield model
void computeDamageOnQuadPoints(GhostType ghost_type);
/// computes strain on quadrature points for phasefield model from
/// displacement gradient from solid mechanics model
void computeStrainOnQuadPoints(GhostType ghost_type);
/// solve the coupled model
void solve(const ID & solid_solver_id = "", const ID & phase_solver_id = "");
private:
/// test the convergence criteria
bool checkConvergence(Array & /*u_new*/, Array & /*u_old*/,
Array & /*d_new*/, Array & /*d_old*/);
protected:
/// callback for the solver, this adds f_{ext} - f_{int} to the residual
void assembleResidual() override;
/// callback for the solver, this adds f_{ext} or f_{int} to the residual
void assembleResidual(const ID & residual_part) override;
bool canSplitResidual() const override { return true; }
/// get the type of matrix needed
MatrixType getMatrixType(const ID & matrix_id) const override;
/// callback for the solver, this assembles different matrices
void assembleMatrix(const ID & matrix_id) override;
/// callback for the solver, this assembles the stiffness matrix
void assembleLumpedMatrix(const ID & matrix_id) override;
/// callback for the model to instantiate the matricess when needed
void initSolver(TimeStepSolverType /*time_step_solver_type*/,
NonLinearSolverType /*non_linear_solver_type*/)
override;
/// callback for the solver, this is called at beginning of solve
void predictor() override;
/// callback for the solver, this is called at end of solve
void corrector() override;
/// 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;
/// solve the coupled model
// void solveStep(const ID & solver_id = "") override;
/// solve a step using a given pre instantiated time step solver and
/// non linear solver with a user defined callback instead of the
/// model itself /!\ This can mess up everything
// void solveStep(SolverCallback & callback, const ID & solver_id = "")
// override;
/* ------------------------------------------------------------------------ */
/* Mass matrix for solid mechanics model */
/* ------------------------------------------------------------------------ */
public:
/// assemble the lumped mass matrix
void assembleMassLumped();
/// assemble the mass matrix for consistent mass resolutions
void assembleMass();
protected:
/// assemble the lumped mass matrix for local and ghost elements
void assembleMassLumped(GhostType ghost_type);
/// assemble the mass matrix for either _ghost or _not_ghost elements
void assembleMass(GhostType ghost_type);
protected:
/* --------------------------------------------------------------------------
*/
TimeStepSolverType getDefaultSolverType() const override;
/* --------------------------------------------------------------------------
*/
ModelSolverOptions
getDefaultSolverOptions(const TimeStepSolverType & type) const override;
public:
bool isDefaultSolverExplicit() { return method == _explicit_lumped_mass; }
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
FEEngine & getFEEngineBoundary(const ID & name = "") override;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// get the solid mechanics model
AKANTU_GET_MACRO(SolidMechanicsModel, *solid, SolidMechanicsModel &);
/// get the contact mechanics model
AKANTU_GET_MACRO(PhaseFieldModel, *phase, PhaseFieldModel &);
/* ------------------------------------------------------------------------ */
/* Dumpable interface */
/* ------------------------------------------------------------------------ */
public:
std::shared_ptr
createNodalFieldReal(const std::string & field_name,
const std::string & group_name,
bool padding_flag) override;
std::shared_ptr
createNodalFieldBool(const std::string & field_name,
const std::string & group_name,
bool padding_flag) override;
std::shared_ptr
createElementalField(const std::string & field_name,
const std::string & group_name, bool padding_flag,
Int spatial_dimension, ElementKind kind) override;
void dump(const std::string & dumper_name) override;
void dump(const std::string & dumper_name, Int step) override;
void dump(const std::string & dumper_name, Real time, Int step) override;
void dump() override;
void dump(Int step) override;
void dump(Real time, Int step) override;
/* ------------------------------------------------------------------------ */
/* Members */
/* ------------------------------------------------------------------------ */
private:
/// solid mechanics model
SolidMechanicsModel * solid{nullptr};
/// phasefield model
PhaseFieldModel * phase{nullptr};
Array * displacement{nullptr};
///
Array * displacement_increment{nullptr};
/// external forces array
Array * external_force{nullptr};
};
} // namespace akantu
#endif /* __AKANTU_COUPLER_SOLID_PHASEFIELD_HH__ */
diff --git a/src/model/phase_field/phase_field_element_filter.hh b/src/model/phase_field/phase_field_element_filter.hh
new file mode 100644
index 000000000..6231cc586
--- /dev/null
+++ b/src/model/phase_field/phase_field_element_filter.hh
@@ -0,0 +1,52 @@
+#include "group_manager.hh"
+#include "phase_field_model.hh"
+
+/* -------------------------------------------------------------------------- */
+namespace akantu {
+
+class PhaseFieldElementFilter : public GroupManager::ClusteringFilter {
+public:
+ PhaseFieldElementFilter(const PhaseFieldModel & model,
+ const Real max_damage = 1.)
+ : model(model), is_unbroken(max_damage) {}
+
+ bool operator()(const Element & el) const override {
+
+ const Array & mat_indexes =
+ model.getPhaseFieldByElement(el.type, el.ghost_type);
+ const Array & mat_loc_num =
+ model.getPhaseFieldLocalNumbering(el.type, el.ghost_type);
+
+ const auto & mat = model.getPhaseField(mat_indexes(el.element));
+
+ Idx el_index = mat_loc_num(el.element);
+ Int nb_quad_per_element =
+ model.getFEEngine("PhaseFieldFEEngine")
+ .getNbIntegrationPoints(el.type, el.ghost_type);
+
+ const Array & damage_array = mat.getDamage(el.type, el.ghost_type);
+
+ AKANTU_DEBUG_ASSERT(nb_quad_per_element * el_index < damage_array.size(),
+ "This quadrature point is out of range");
+
+ const Real * element_damage =
+ damage_array.data() + nb_quad_per_element * el_index;
+
+ UInt unbroken_quads = std::count_if(
+ element_damage, element_damage + nb_quad_per_element, is_unbroken);
+
+ return (unbroken_quads > 0);
+ }
+
+private:
+ struct IsUnbrokenFunctor {
+ IsUnbrokenFunctor(const Real & max_damage) : max_damage(max_damage) {}
+ bool operator()(const Real & x) const { return x > max_damage; }
+ const Real max_damage;
+ };
+
+ const PhaseFieldModel & model;
+ const IsUnbrokenFunctor is_unbroken;
+};
+
+}
diff --git a/src/model/phase_field/phase_field_model.cc b/src/model/phase_field/phase_field_model.cc
index 2581fdf52..03be931d5 100644
--- a/src/model/phase_field/phase_field_model.cc
+++ b/src/model/phase_field/phase_field_model.cc
@@ -1,583 +1,744 @@
/**
* 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 .
*/
/* -------------------------------------------------------------------------- */
#include "phase_field_model.hh"
+#include "aka_common.hh"
#include "dumpable_inline_impl.hh"
#include "element_synchronizer.hh"
#include "fe_engine_template.hh"
#include "generalized_trapezoidal.hh"
#include "group_manager_inline_impl.hh"
#include "integrator_gauss.hh"
#include "mesh.hh"
#include "parser.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
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
PhaseFieldModel::PhaseFieldModel(Mesh & mesh, Int dim, const ID & id,
+ std::shared_ptr dof_manager,
ModelType model_type)
: Model(mesh, model_type, dim, id),
phasefield_index("phasefield index", id),
phasefield_local_numbering("phasefield local numbering", id) {
AKANTU_DEBUG_IN();
+ this->initDOFManager(std::move(dof_manager));
+
this->registerFEEngineObject("PhaseFieldFEEngine", mesh,
Model::spatial_dimension);
this->mesh.registerDumper("phase_field", id, true);
this->mesh.addDumpMesh(mesh, Model::spatial_dimension, _not_ghost,
_ek_regular);
phasefield_selector =
std::make_shared(phasefield_index);
- this->initDOFManager();
-
this->registerDataAccessor(*this);
if (this->mesh.isDistributed()) {
auto & synchronizer = this->mesh.getElementSynchronizer();
+ this->registerSynchronizer(synchronizer,
+ SynchronizationTag::_phasefield_id);
this->registerSynchronizer(synchronizer, SynchronizationTag::_pfm_damage);
- this->registerSynchronizer(synchronizer, SynchronizationTag::_pfm_driving);
- this->registerSynchronizer(synchronizer, SynchronizationTag::_pfm_history);
- this->registerSynchronizer(synchronizer, SynchronizationTag::_pfm_energy);
+ this->registerSynchronizer(synchronizer, SynchronizationTag::_for_dump);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
PhaseFieldModel::~PhaseFieldModel() = default;
/* -------------------------------------------------------------------------- */
MatrixType PhaseFieldModel::getMatrixType(const ID & matrix_id) const {
if (matrix_id == "K") {
return _symmetric;
}
return _mt_not_defined;
}
/* -------------------------------------------------------------------------- */
void PhaseFieldModel::initModel() {
auto & fem = this->getFEEngine();
fem.initShapeFunctions(_not_ghost);
fem.initShapeFunctions(_ghost);
}
/* -------------------------------------------------------------------------- */
void PhaseFieldModel::initFullImpl(const ModelOptions & options) {
phasefield_index.initialize(mesh, _element_kind = _ek_not_defined,
- _default_value = UInt(-1),
+ _default_value = Idx(-1),
_with_nb_element = true);
phasefield_local_numbering.initialize(mesh, _element_kind = _ek_not_defined,
_with_nb_element = true);
Model::initFullImpl(options);
// initialize the phasefields
if (!this->parser.getLastParsedFile().empty()) {
this->instantiatePhaseFields();
this->initPhaseFields();
}
this->initBC(*this, *damage, *external_force);
}
/* -------------------------------------------------------------------------- */
PhaseField &
PhaseFieldModel::registerNewPhaseField(const ParserSection & section) {
std::string phase_name;
std::string phase_type = section.getName();
std::string opt_param = section.getOption();
try {
std::string tmp = section.getParameter("name");
phase_name = tmp; /** this can seam weird, but there is an ambiguous
* operator overload that i couldn't solve. @todo remove
* the weirdness of this code
*/
} catch (debug::Exception &) {
AKANTU_ERROR("A phasefield of type \'"
<< phase_type
<< "\' in the input file has been defined without a name!");
}
PhaseField & phase =
this->registerNewPhaseField(phase_name, phase_type, opt_param);
phase.parseSection(section);
return phase;
}
/* -------------------------------------------------------------------------- */
PhaseField & PhaseFieldModel::registerNewPhaseField(const ID & phase_name,
const ID & phase_type,
const ID & opt_param) {
AKANTU_DEBUG_ASSERT(phasefields_names_to_id.find(phase_name) ==
phasefields_names_to_id.end(),
"A phasefield with this name '"
<< phase_name << "' has already been registered. "
<< "Please use unique names for phasefields");
- UInt phase_count = phasefields.size();
+ Int phase_count = phasefields.size();
phasefields_names_to_id[phase_name] = phase_count;
std::stringstream sstr_phase;
sstr_phase << this->id << ":" << phase_count << ":" << phase_type;
ID mat_id = sstr_phase.str();
std::unique_ptr phase = PhaseFieldFactory::getInstance().allocate(
phase_type, opt_param, *this, mat_id);
phasefields.push_back(std::move(phase));
return *(phasefields.back());
}
/* -------------------------------------------------------------------------- */
void PhaseFieldModel::instantiatePhaseFields() {
ParserSection model_section;
bool is_empty;
std::tie(model_section, is_empty) = this->getParserSection();
if (not is_empty) {
auto model_phasefields =
model_section.getSubSections(ParserType::_phasefield);
for (const auto & section : model_phasefields) {
this->registerNewPhaseField(section);
}
}
auto sub_sections = this->parser.getSubSections(ParserType::_phasefield);
for (const auto & section : sub_sections) {
this->registerNewPhaseField(section);
}
if (phasefields.empty()) {
AKANTU_EXCEPTION("No phasefields where instantiated for the model"
<< getID());
}
are_phasefields_instantiated = true;
}
/* -------------------------------------------------------------------------- */
void PhaseFieldModel::initPhaseFields() {
AKANTU_DEBUG_ASSERT(phasefields.size() != 0, "No phasefield to initialize !");
if (!are_phasefields_instantiated) {
instantiatePhaseFields();
}
this->assignPhaseFieldToElements();
for (auto & phasefield : phasefields) {
/// init internals properties
phasefield->initPhaseField();
}
-
- this->synchronize(SynchronizationTag::_smm_init_mat);
}
/* -------------------------------------------------------------------------- */
void PhaseFieldModel::assignPhaseFieldToElements(
const ElementTypeMapArray * filter) {
for_each_element(
mesh,
[&](auto && element) {
Int phase_index = (*phasefield_selector)(element);
AKANTU_DEBUG_ASSERT(
phase_index < Int(phasefields.size()),
"The phasefield selector returned an index that does not exists");
phasefield_index(element) = phase_index;
},
_element_filter = filter, _ghost_type = _not_ghost);
for_each_element(
mesh,
[&](auto && element) {
auto phase_index = phasefield_index(element);
auto index = phasefields[phase_index]->addElement(element);
phasefield_local_numbering(element) = index;
},
_element_filter = filter, _ghost_type = _not_ghost);
// synchronize the element phasefield arrays
- this->synchronize(SynchronizationTag::_material_id);
+ this->synchronize(SynchronizationTag::_phasefield_id);
}
/* -------------------------------------------------------------------------- */
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(getFEEngineClassBoundary(name));
}
+/* -------------------------------------------------------------------------- */
+TimeStepSolverType PhaseFieldModel::getDefaultSolverType() const {
+ return TimeStepSolverType::_static;
+}
+
/* -------------------------------------------------------------------------- */
std::tuple
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::_linear;
+ 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() {
+ AKANTU_DEBUG_IN();
+
+ Real energy = 0.;
+ for (auto & phasefield : phasefields) {
+ energy += phasefield->getEnergy();
+ }
+
+ /// reduction sum over all processors
+ mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum);
+
+ AKANTU_DEBUG_OUT();
+ return energy;
+}
+
+/* -------------------------------------------------------------------------- */
+Real PhaseFieldModel::getEnergy(ElementType type, Idx index) {
+ AKANTU_DEBUG_IN();
+
+ Idx phase_index = this->phasefield_index(type, _not_ghost)(index);
+ Idx phase_loc_num = this->phasefield_local_numbering(type, _not_ghost)(index);
+ Real energy = this->phasefields[phase_index]->getEnergy(
+ Element{type, phase_loc_num, _not_ghost});
+
+ AKANTU_DEBUG_OUT();
+ return energy;
+}
+
+/* -------------------------------------------------------------------------- */
+Real PhaseFieldModel::getEnergy(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 += getEnergy(el);
+ }
+ }
+
+ /// reduction sum over all processors
+ mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum);
+
+ return energy;
+}
+
/* -------------------------------------------------------------------------- */
void PhaseFieldModel::beforeSolveStep() {
for (auto & phasefield : phasefields) {
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);
- dam -= prev_dam;
prev_dam = dam;
}
}
/* -------------------------------------------------------------------------- */
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 (auto & phasefield : phasefields) {
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();
- // communicate the driving forces
- AKANTU_DEBUG_INFO("Send data for residual assembly");
- this->asynchronousSynchronize(SynchronizationTag::_pfm_driving);
+ this->synchronize(SynchronizationTag::_pfm_damage);
+
+ for (auto & phasefield : phasefields) {
+ phasefield->computeAllDrivingForces(_not_ghost);
+ }
// assemble the forces due to local driving forces
AKANTU_DEBUG_INFO("Assemble residual for local elements");
for (auto & phasefield : phasefields) {
phasefield->assembleInternalForces(_not_ghost);
}
- // finalize communications
- AKANTU_DEBUG_INFO("Wait distant driving forces");
- this->waitEndSynchronize(SynchronizationTag::_pfm_driving);
-
- // assemble the residual due to ghost elements
+ // assemble the forces due to local driving forces
AKANTU_DEBUG_INFO("Assemble residual for ghost elements");
+ for (auto & phasefield : phasefields) {
+ phasefield->assembleInternalForces(_ghost);
+ }
}
/* -------------------------------------------------------------------------- */
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 & 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::_pfm_damage: {
- size += nb_nodes_per_element * sizeof(Real); // damage
- break;
- }
- case SynchronizationTag::_pfm_driving: {
- size += getNbIntegrationPoints(elements) * sizeof(Real);
+ case SynchronizationTag::_phasefield_id: {
+ size += elements.size() * sizeof(Int);
break;
}
- case SynchronizationTag::_pfm_history: {
- size += getNbIntegrationPoints(elements) * sizeof(Real);
+ case SynchronizationTag::_for_dump: {
+ // damage
+ size += nb_nodes_per_element * sizeof(Real);
break;
}
- case SynchronizationTag::_pfm_energy: {
- size += getNbIntegrationPoints(elements) * sizeof(Real);
+ case SynchronizationTag::_pfm_damage: {
+ size += nb_nodes_per_element * sizeof(Real);
break;
}
default: {
AKANTU_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
return size;
}
/* -------------------------------------------------------------------------- */
-void PhaseFieldModel::packData(CommunicationBuffer & /*buffer*/,
- const Array & /*elements*/,
- const SynchronizationTag & /*tag*/) const {}
+void PhaseFieldModel::packData(CommunicationBuffer & buffer,
+ const Array & elements,
+ const SynchronizationTag & tag) const {
+ switch (tag) {
+ case SynchronizationTag::_phasefield_id: {
+ packElementalDataHelper(phasefield_index, buffer, elements, false,
+ getFEEngine());
+ break;
+ }
+ 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);
+ }
+ }
+}
/* -------------------------------------------------------------------------- */
-void PhaseFieldModel::unpackData(CommunicationBuffer & /*buffer*/,
- const Array & /*elements*/,
- const SynchronizationTag & /*tag*/) {}
+void PhaseFieldModel::unpackData(CommunicationBuffer & buffer,
+ const Array & elements,
+ const SynchronizationTag & tag) {
+ AKANTU_DEBUG_IN();
+
+ switch (tag) {
+ case SynchronizationTag::_phasefield_id: {
+ for (auto && element : elements) {
+ Idx recv_phase_index;
+ buffer >> recv_phase_index;
+ Idx & phase_index = phasefield_index(element);
+ if (phase_index != Idx(-1)) {
+ continue;
+ }
+
+ // add ghosts element to the correct phasefield
+ phase_index = recv_phase_index;
+ Idx index = phasefields[phase_index]->addElement(element);
+ phasefield_local_numbering(element) = index;
+ }
+ break;
+ }
+ 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);
+ }
+ }
+
+ AKANTU_DEBUG_OUT();
+}
/* -------------------------------------------------------------------------- */
Int PhaseFieldModel::getNbData(const Array & indexes,
const SynchronizationTag & tag) const {
- UInt size = 0;
- UInt nb_nodes = indexes.size();
+ 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 & indexes,
const SynchronizationTag & tag) const {
- for (auto index : indexes) {
- switch (tag) {
- case SynchronizationTag::_pfm_damage: {
- buffer << (*damage)(index);
- break;
- }
- default: {
- AKANTU_ERROR("Unknown ghost synchronization tag : " << tag);
- }
- }
+ 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 & indexes,
const SynchronizationTag & tag) {
- for (auto index : indexes) {
- switch (tag) {
- case SynchronizationTag::_pfm_damage: {
- buffer >> (*damage)(index);
- break;
- }
- default: {
- AKANTU_ERROR("Unknown ghost synchronization tag : " << 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
PhaseFieldModel::createNodalFieldBool(const std::string & field_name,
const std::string & group_name,
bool /*unused*/) {
std::map *> uint_nodal_fields;
uint_nodal_fields["blocked_dofs"] = blocked_dofs.get();
return mesh.createNodalField(uint_nodal_fields[field_name], group_name);
std::shared_ptr field;
return field;
}
/* -------------------------------------------------------------------------- */
std::shared_ptr
PhaseFieldModel::createNodalFieldReal(const std::string & field_name,
const std::string & group_name,
bool /*unused*/) {
std::map *> 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 field;
return field;
}
/* -------------------------------------------------------------------------- */
std::shared_ptr 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(
mesh.getConnectivities(), group_name, this->spatial_dimension,
element_kind);
}
std::shared_ptr field;
return field;
}
+/* -------------------------------------------------------------------------- */
+ElementTypeMapArray &
+PhaseFieldModel::flattenInternal(const std::string & field_name,
+ ElementKind kind, const GhostType ghost_type) {
+ auto key = std::make_pair(field_name, kind);
+
+ ElementTypeMapArray * internal_flat;
+
+ auto it = this->registered_internals.find(key);
+ if (it == this->registered_internals.end()) {
+ auto internal =
+ std::make_unique>(field_name, this->id);
+
+ internal_flat = internal.get();
+ this->registered_internals[key] = std::move(internal);
+ } else {
+ internal_flat = it->second.get();
+ }
+
+ for (auto type :
+ mesh.elementTypes(Model::spatial_dimension, ghost_type, kind)) {
+ if (internal_flat->exists(type, ghost_type)) {
+ auto & internal = (*internal_flat)(type, ghost_type);
+ internal.resize(0);
+ }
+ }
+
+ for (auto & phasefield : phasefields) {
+ if (phasefield->isInternal(field_name, kind)) {
+ phasefield->flattenInternal(field_name, *internal_flat, ghost_type, kind);
+ }
+ }
+
+ return *internal_flat;
+}
+
+/* -------------------------------------------------------------------------- */
+void PhaseFieldModel::inflateInternal(const std::string & field_name,
+ const ElementTypeMapArray & field,
+ ElementKind kind, GhostType ghost_type) {
+
+ for (auto & phasefield : phasefields) {
+ if (phasefield->isInternal(field_name, kind)) {
+ phasefield->inflateInternal(field_name, field, ghost_type, kind);
+ } else {
+ AKANTU_ERROR("A internal of name \'"
+ << field_name
+ << "\' has not been defined in the phasefield");
+ }
+ }
+}
+
/* -------------------------------------------------------------------------- */
void PhaseFieldModel::printself(std::ostream & stream, int indent) const {
std::string space(indent, AKANTU_INDENT);
stream << space << "Phase Field Model [" << std::endl;
stream << space << " + id : " << id << std::endl;
stream << space << " + spatial dimension : " << Model::spatial_dimension
<< std::endl;
stream << space << " + fem [" << std::endl;
getFEEngine().printself(stream, indent + 2);
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << " + nodals information [" << std::endl;
damage->printself(stream, indent + 2);
external_force->printself(stream, indent + 2);
internal_force->printself(stream, indent + 2);
blocked_dofs->printself(stream, indent + 2);
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << " + phasefield information [" << std::endl;
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << "]" << std::endl;
}
} // namespace akantu
diff --git a/src/model/phase_field/phase_field_model.hh b/src/model/phase_field/phase_field_model.hh
index 2e36c3789..7f178fa32 100644
--- a/src/model/phase_field/phase_field_model.hh
+++ b/src/model/phase_field/phase_field_model.hh
@@ -1,327 +1,367 @@
/**
* 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 .
*/
/* -------------------------------------------------------------------------- */
#include "boundary_condition.hh"
#include "data_accessor.hh"
#include "fe_engine.hh"
#include "model.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_PHASE_FIELD_MODEL_HH__
#define __AKANTU_PHASE_FIELD_MODEL_HH__
namespace akantu {
class PhaseField;
class PhaseFieldSelector;
template class IntegratorGauss;
template class ShapeLagrange;
} // namespace akantu
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
class PhaseFieldModel : public Model,
public DataAccessor,
public DataAccessor,
public BoundaryCondition {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
using FEEngineType = FEEngineTemplate;
PhaseFieldModel(Mesh & mesh, Int dim = _all_dimensions,
const ID & id = "phase_field_model",
+ std::shared_ptr dof_manager = nullptr,
ModelType model_type = ModelType::_phase_field_model);
~PhaseFieldModel() override;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
protected:
/// generic function to initialize everything ready for explicit dynamics
void initFullImpl(const ModelOptions & options) override;
/// initialize all internal array for phasefields
void initPhaseFields();
/// allocate all vectors
void initSolver(TimeStepSolverType /*unused*/,
NonLinearSolverType /*unused*/) override;
/// initialize the model
void initModel() override;
/// predictor
void predictor() override;
/// corrector
void corrector() override;
/// compute the heat flux
void assembleResidual() override;
/// get the type of matrix needed
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;
+ /// function to print the containt of the class
+ void printself(std::ostream & stream, int indent = 0) const override;
+
+protected:
+ /* ------------------------------------------------------------------------ */
+ TimeStepSolverType getDefaultSolverType() const override;
+
std::tuple
getDefaultSolverID(const AnalysisMethod & method) override;
ModelSolverOptions
getDefaultSolverOptions(const TimeStepSolverType & type) const override;
- /// function to print the containt of the class
- void printself(std::ostream & stream, int indent = 0) const override;
-
/* ------------------------------------------------------------------------ */
/* Materials (phase_field_model.cc) */
/* ------------------------------------------------------------------------ */
public:
/// register an empty phasefield of a given type
PhaseField & registerNewPhaseField(const ID & phase_name,
const ID & phase_type,
const ID & opt_param);
/// reassigns phasefields depending on the phasefield selector
void reassignPhaseField();
protected:
/// register a phasefield in the dynamic database
PhaseField & registerNewPhaseField(const ParserSection & phase_section);
/// read the phasefield files to instantiate all the phasefields
void instantiatePhaseFields();
/// set the element_id_by_phasefield and add the elements to the good
/// phasefields
void
assignPhaseFieldToElements(const ElementTypeMapArray * 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);
/* ------------------------------------------------------------------------ */
/* 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;
/* ------------------------------------------------------------------------ */
/* Data Accessor inherited members */
/* ------------------------------------------------------------------------ */
public:
Int getNbData(const Array & elements,
const SynchronizationTag & tag) const override;
void packData(CommunicationBuffer & buffer, const Array & elements,
const SynchronizationTag & tag) const override;
void unpackData(CommunicationBuffer & buffer, const Array & elements,
const SynchronizationTag & tag) override;
Int getNbData(const Array & indexes,
const SynchronizationTag & tag) const override;
void packData(CommunicationBuffer & buffer, const Array & indexes,
const SynchronizationTag & tag) const override;
void unpackData(CommunicationBuffer & buffer, const Array & indexes,
const SynchronizationTag & tag) override;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
+protected:
+ /// Compute dissipated energy for an element type and phasefield index
+ Real getEnergy(ElementType type, Idx index);
+
public:
/// return the damage array
AKANTU_GET_MACRO_DEREF_PTR(Damage, damage);
AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Damage, 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::force vector (external forces)
Array & getForce() {
AKANTU_DEBUG_WARNING("getForce was maintained for backward compatibility, "
"use getExternalForce instead");
return *external_force;
}
/// get the PhaseFieldModel::blocked_dofs vector
AKANTU_GET_MACRO_DEREF_PTR(BlockedDOFs, blocked_dofs);
/// get an iterable on the phasefields
inline decltype(auto) getPhaseFields();
/// get an iterable on the phasefields
inline decltype(auto) getPhaseFields() const;
/// get a particular phasefield (by phasefield index)
- inline PhaseField & getPhaseField(UInt mat_index);
+ inline PhaseField & getPhaseField(Idx mat_index);
/// get a particular phasefield (by phasefield index)
- inline const PhaseField & getPhaseField(UInt mat_index) const;
+ inline const PhaseField & getPhaseField(Idx mat_index) const;
/// get a particular phasefield (by phasefield name)
inline PhaseField & getPhaseField(const std::string & name);
/// get a particular phasefield (by phasefield name)
inline const PhaseField & getPhaseField(const std::string & name) const;
/// get a particular phasefield id from is name
- inline Int getPhaseFieldIndex(const std::string & name) const;
+ inline Idx getPhaseFieldIndex(const std::string & name) const;
/// give the number of phasefields
- inline Int getNbPhaseFields() const { return phasefields.size(); }
+ inline Idx getNbPhaseFields() const { return phasefields.size(); }
/// give the phasefield internal index from its id
- Int getInternalIndexFromID(const ID & id) const;
+ Idx getInternalIndexFromID(const ID & id) const;
+
+ /**
+ * @brief Returns the total dissipated energy
+ *
+ */
+ Real getEnergy();
+
+ /// Compute dissipated energy for an individual element
+ Real getEnergy(const Element & element) {
+ return getEnergy(element.type, element.element);
+ }
+
+ /// Compute dissipated energy for an element group
+ Real getEnergy(const ID & group_id);
AKANTU_GET_MACRO(PhaseFieldByElement, phasefield_index,
const ElementTypeMapArray &);
AKANTU_GET_MACRO(PhaseFieldLocalNumbering, phasefield_local_numbering,
const ElementTypeMapArray &);
/// vectors containing local material element index for each global element
/// index
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(PhaseFieldByElement, phasefield_index,
Idx);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE(PhaseFieldByElement, phasefield_index, Idx);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(PhaseFieldLocalNumbering,
phasefield_local_numbering, Idx);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE(PhaseFieldLocalNumbering,
phasefield_local_numbering, Idx);
AKANTU_GET_MACRO_NOT_CONST(PhaseFieldSelector, *phasefield_selector,
PhaseFieldSelector &);
AKANTU_SET_MACRO(PhaseFieldSelector, phasefield_selector,
std::shared_ptr);
FEEngine & getFEEngineBoundary(const ID & name = "") override;
/* ------------------------------------------------------------------------ */
/* Dumpable Interface */
/* ------------------------------------------------------------------------ */
public:
std::shared_ptr
createNodalFieldReal(const std::string & field_name,
const std::string & group_name,
bool padding_flag) override;
std::shared_ptr
createNodalFieldBool(const std::string & field_name,
const std::string & group_name,
bool padding_flag) override;
std::shared_ptr
createElementalField(const std::string & field_name,
const std::string & group_name, bool padding_flag,
Int spatial_dimension, ElementKind kind) override;
+ //! flatten a given phasefield internal field
+ ElementTypeMapArray &
+ flattenInternal(const std::string & field_name, ElementKind kind,
+ GhostType ghost_type = _not_ghost);
+
+ //! inverse operation of the flatten
+ void inflateInternal(const std::string & field_name,
+ const ElementTypeMapArray & field,
+ ElementKind kind, GhostType ghost_type = _not_ghost);
+
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
/// number of iterations
Int n_iter;
/// damage array
std::unique_ptr> damage;
/// damage array at the previous time step
std::unique_ptr> previous_damage;
/// boundary vector
std::unique_ptr> blocked_dofs;
/// external force vector
std::unique_ptr> external_force;
/// residuals array
std::unique_ptr> internal_force;
/// Arrays containing the phasefield index for each element
ElementTypeMapArray phasefield_index;
/// Arrays containing the position in the element filter of the phasefield
/// (phasefield's local numbering)
ElementTypeMapArray phasefield_local_numbering;
/// class defining of to choose a phasefield
std::shared_ptr phasefield_selector;
/// mapping between phasefield name and phasefield internal id
std::map phasefields_names_to_id;
/// list of used phasefields
std::vector> phasefields;
+ using flatten_internal_map =
+ std::map,
+ std::unique_ptr>>;
+
+ /// tells if the phasefields are instantiated
+ flatten_internal_map registered_internals;
+
/// tells if the phasefield are instantiated
bool are_phasefields_instantiated{false};
};
} // namespace akantu
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#include "parser.hh"
#include "phasefield.hh"
-#include "phase_field_model_inline_impl.cc"
+#include "phase_field_model_inline_impl.hh"
/* -------------------------------------------------------------------------- */
#endif
diff --git a/src/model/phase_field/phase_field_model_inline_impl.cc b/src/model/phase_field/phase_field_model_inline_impl.hh
similarity index 96%
rename from src/model/phase_field/phase_field_model_inline_impl.cc
rename to src/model/phase_field/phase_field_model_inline_impl.hh
index 431a1319e..afc009270 100644
--- a/src/model/phase_field/phase_field_model_inline_impl.cc
+++ b/src/model/phase_field/phase_field_model_inline_impl.hh
@@ -1,90 +1,90 @@
/**
* 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_named_argument.hh"
#include "phasefield_selector.hh"
#include "phasefield_selector_tmpl.hh"
#include "solid_mechanics_model.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_PHASE_FIELD_MODEL_INLINE_IMPL_CC__
#define __AKANTU_PHASE_FIELD_MODEL_INLINE_IMPL_CC__
namespace akantu {
/* -------------------------------------------------------------------------- */
inline decltype(auto) PhaseFieldModel::getPhaseFields() {
return make_dereference_adaptor(phasefields);
}
/* -------------------------------------------------------------------------- */
inline decltype(auto) PhaseFieldModel::getPhaseFields() const {
return make_dereference_adaptor(phasefields);
}
/* -------------------------------------------------------------------------- */
-inline PhaseField & PhaseFieldModel::getPhaseField(UInt mat_index) {
+inline PhaseField & PhaseFieldModel::getPhaseField(Idx mat_index) {
AKANTU_DEBUG_ASSERT(mat_index < phasefields.size(),
"The model " << id << " has no phasefield no "
<< mat_index);
return *phasefields[mat_index];
}
/* -------------------------------------------------------------------------- */
-inline const PhaseField & PhaseFieldModel::getPhaseField(UInt mat_index) const {
+inline const PhaseField & PhaseFieldModel::getPhaseField(Idx mat_index) const {
AKANTU_DEBUG_ASSERT(mat_index < phasefields.size(),
"The model " << id << " has no phasefield no "
<< mat_index);
return *phasefields[mat_index];
}
/* -------------------------------------------------------------------------- */
inline PhaseField & PhaseFieldModel::getPhaseField(const std::string & name) {
auto it = phasefields_names_to_id.find(name);
AKANTU_DEBUG_ASSERT(it != phasefields_names_to_id.end(),
"The model " << id << " has no phasefield named "
<< name);
return *phasefields[it->second];
}
/* -------------------------------------------------------------------------- */
inline Int PhaseFieldModel::getPhaseFieldIndex(const std::string & name) const {
auto it = phasefields_names_to_id.find(name);
AKANTU_DEBUG_ASSERT(it != phasefields_names_to_id.end(),
"The model " << id << " has no phasefield named "
<< name);
return it->second;
}
/* -------------------------------------------------------------------------- */
inline const PhaseField &
PhaseFieldModel::getPhaseField(const std::string & name) const {
auto it = phasefields_names_to_id.find(name);
AKANTU_DEBUG_ASSERT(it != phasefields_names_to_id.end(),
"The model " << id << " has no phasefield named "
<< name);
return *phasefields[it->second];
}
/* -------------------------------------------------------------------------- */
} // namespace akantu
#endif /* __AKANTU_PHASE_FIELD_MODEL_INLINE_IMPL_CC__ */
diff --git a/src/model/phase_field/phasefield.cc b/src/model/phase_field/phasefield.cc
index 3bff8a6c9..94ff1f5fd 100644
--- a/src/model/phase_field/phasefield.cc
+++ b/src/model/phase_field/phasefield.cc
@@ -1,274 +1,409 @@
/**
* 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("damage", *this, dim, fe_engine, this->element_filter),
+ 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) {
+ 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(PhaseFieldModel & model, const ID & id)
- : PhaseField(model, model.getSpatialDimension(), model.getMesh(),
- model.getFEEngine(), id) {}
-
/* -------------------------------------------------------------------------- */
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");
- damage.initialize(1);
+ 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(type, ghost_type), nt_driving_force, type,
- ghost_type, elem_filter);
+ fem.computeNtb(driving_force_vect, nt_driving_force, type, ghost_type,
+ elem_filter);
- Array int_nt_driving_force(nb_element, nb_nodes_per_element);
+ 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);
+ 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 cfaf66495..85240799f 100644
--- a/src/model/phase_field/phasefield.hh
+++ b/src/model/phase_field/phasefield.hh
@@ -1,298 +1,357 @@
/**
* 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 UInt addElement(ElementType type, UInt element, GhostType ghost_type);
- inline UInt addElement(const Element & element);
+ 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, Real);
+ AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Damage, damage_on_qpoints, Real);
- AKANTU_GET_MACRO_NOT_CONST(Damage, damage, ElementTypeMapArray &);
- AKANTU_GET_MACRO(Damage, damage, const ElementTypeMapArray &);
+ 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, ElementKind element_kind) const;
+ 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_int;
std::map