diff --git a/examples/phase_field/phase_field_notch.cc b/examples/phase_field/phase_field_notch.cc index 2b114e21f..8b1850821 100644 --- a/examples/phase_field/phase_field_notch.cc +++ b/examples/phase_field/phase_field_notch.cc @@ -1,182 +1,182 @@ /** * @file phase_field_notch.cc * * @author Mohit Pundir * * @date creation: Tue Oct 02 2018 * @date last modification: Wed Apr 07 2021 * * @brief Example of phase field model * * * @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 "coupler_solid_phasefield.hh" #include "non_linear_solver.hh" #include "phase_field_model.hh" #include "solid_mechanics_model.hh" #include "group_manager.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 UInt spatial_dimension = 2; /* -------------------------------------------------------------------------- */ 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)); UInt el_index = mat_loc_num(el.element); UInt 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.storage() + 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; }; 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(); UInt nbSteps = 1500; Real increment = 1e-5; auto start_time = clk::now(); for (UInt s = 1; s < nbSteps; ++s) { if (s >= 500) { increment = 1.e-6; } 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(); auto energy = phase.getEnergy(); if (s % 100 == 0) { model.dump(); } } - Real damage_limit = 0.08; - auto global_nb_clusters = - mesh.createClusters(spatial_dimension, "crack", PhaseFieldElementFilter(phase, damage_limit)); + // Real damage_limit = 0.08; + // auto global_nb_clusters = + // mesh.createClusters(spatial_dimension, "crack", PhaseFieldElementFilter(phase, damage_limit)); - - auto nb_fragment = mesh.getNbElementGroups(spatial_dimension); + // + // auto nb_fragment = mesh.getNbElementGroups(spatial_dimension); - model.dumpGroup("crack_0"); - - std::cout << global_nb_clusters << std::endl; - std::cout << nb_fragment << std::endl; + // model.dumpGroup("crack"); + // + // std::cout << global_nb_clusters << std::endl; + // std::cout << nb_fragment << std::endl; finalize(); return EXIT_SUCCESS; } diff --git a/examples/phase_field/phase_field_parallel.cc b/examples/phase_field/phase_field_parallel.cc index 3cae1eb89..cee28d2a7 100644 --- a/examples/phase_field/phase_field_parallel.cc +++ b/examples/phase_field/phase_field_parallel.cc @@ -1,133 +1,133 @@ /** * @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 * * * @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 "non_linear_solver.hh" #include "phase_field_model.hh" #include "solid_mechanics_model.hh" #include "group_manager.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 UInt spatial_dimension = 2; int main(int argc, char * argv[]) { initialize("material_notch.dat", argc, argv); // create mesh Mesh mesh(spatial_dimension); const auto & comm = Communicator::getStaticCommunicator(); Int prank = comm.whoAmI(); if (prank == 0) { // Read the mesh mesh.read("square_notch.msh"); } mesh.distribute(); CouplerSolidPhaseField coupler(mesh); auto & model = coupler.getSolidMechanicsModel(); auto & phase = coupler.getPhaseFieldModel(); model.initFull(_analysis_method = _static); auto && mat_selector = std::make_shared>("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"); - phase.setBaseName("phase_notch_parallel"); - //model.addDumpField("stress"); - //model.addDumpField("grad_u"); - //model.addDumpFieldVector("displacement"); - phase.addDumpField("damage"); + model.setBaseName("phase_notch_parallel"); + model.addDumpField("stress"); + model.addDumpField("grad_u"); + model.addDumpFieldVector("displacement"); + model.addDumpField("damage"); if (mesh.isDistributed()) { //phase.addDumpField("partitions"); } phase.dump(); UInt nbSteps = 1500; Real increment = 1e-5; auto start_time = clk::now(); for (UInt s = 1; s < nbSteps; ++s) { if (s >= 500) { increment = 1.e-6; } 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(); if (s % 100 == 0) { - phase.dump(); + model.dump(); } } finalize(); return EXIT_SUCCESS; } diff --git a/python/py_dof_manager.cc b/python/py_dof_manager.cc index 1a073da25..e268e7586 100644 --- a/python/py_dof_manager.cc +++ b/python/py_dof_manager.cc @@ -1,216 +1,216 @@ /* * 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 "py_dof_manager.hh" #include "py_aka_array.hh" #include "py_akantu_pybind11_compatibility.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { namespace { class PySolverCallback : public SolverCallback { public: using SolverCallback::SolverCallback; /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) const override { // NOLINTNEXTLINE PYBIND11_OVERRIDE_PURE(MatrixType, SolverCallback, getMatrixType, matrix_id); } /// callback to assemble a Matrix void assembleMatrix(const ID & matrix_id) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE_PURE(void, SolverCallback, assembleMatrix, matrix_id); } /// callback to assemble a lumped Matrix void assembleLumpedMatrix(const ID & matrix_id) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE_PURE(void, SolverCallback, assembleLumpedMatrix, matrix_id); } /// callback to assemble the residual (rhs) void assembleResidual() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE_PURE(void, SolverCallback, assembleResidual); } /// callback for the predictor (in case of dynamic simulation) void predictor() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, SolverCallback, predictor); } /// callback for the corrector (in case of dynamic simulation) void corrector() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, SolverCallback, corrector); } void beforeSolveStep() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, SolverCallback, beforeSolveStep); } void afterSolveStep(bool converged) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, SolverCallback, afterSolveStep, converged); } }; class PyInterceptSolverCallback : public InterceptSolverCallback { public: using InterceptSolverCallback::InterceptSolverCallback; MatrixType getMatrixType(const ID & matrix_id) const override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(MatrixType, InterceptSolverCallback, getMatrixType, matrix_id); } void assembleMatrix(const ID & matrix_id) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, InterceptSolverCallback, assembleMatrix, matrix_id); } /// callback to assemble a lumped Matrix void assembleLumpedMatrix(const ID & matrix_id) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, InterceptSolverCallback, assembleLumpedMatrix, matrix_id); } void assembleResidual() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, InterceptSolverCallback, assembleResidual); } void predictor() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, InterceptSolverCallback, predictor); } void corrector() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, InterceptSolverCallback, corrector); } void beforeSolveStep() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, InterceptSolverCallback, beforeSolveStep); } void afterSolveStep(bool converged) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, InterceptSolverCallback, afterSolveStep, converged); } }; } // namespace /* -------------------------------------------------------------------------- */ void register_dof_manager(py::module & mod) { py::class_>(mod, "DOFManager") .def("getMatrix", &DOFManager::getMatrix, py::return_value_policy::reference) .def( "getNewMatrix", [](DOFManager & self, const std::string & name, const std::string & matrix_to_copy_id) -> decltype(auto) { return self.getNewMatrix(name, matrix_to_copy_id); }, py::return_value_policy::reference) .def( "getResidual", [](DOFManager & self) -> decltype(auto) { return self.getResidual(); }, py::return_value_policy::reference) .def("getArrayPerDOFs", &DOFManager::getArrayPerDOFs) .def( "hasMatrix", [](DOFManager & self, const ID & name) -> bool { return self.hasMatrix(name); }, py::arg("name")) .def("assembleToResidual", &DOFManager::assembleToResidual, py::arg("dof_id"), py::arg("array_to_assemble"), py::arg("scale_factor") = 1.) .def("assembleToLumpedMatrix", &DOFManager::assembleToLumpedMatrix, py::arg("dof_id"), py::arg("array_to_assemble"), py::arg("lumped_mtx"), py::arg("scale_factor") = 1.) .def("assemblePreassembledMatrix", &DOFManager::assemblePreassembledMatrix, py::arg("matrix_id"), py::arg("terms")) .def("zeroResidual", &DOFManager::zeroResidual); - py::class_(mod, "NonLinearSolver") + py::class_(mod, "NonLinearSolver") .def( "set", [](NonLinearSolver & self, const std::string & id, const Real & val) { if (id == "max_iterations") { self.set(id, int(val)); } else { self.set(id, val); } }) .def("set", [](NonLinearSolver & self, const std::string & id, const SolveConvergenceCriteria & val) { self.set(id, val); }); py::class_(mod, "TimeStepSolver") .def("getIntegrationScheme", &TimeStepSolver::getIntegrationScheme); py::class_(mod, "SolverCallback") .def(py::init_alias()) .def("getMatrixType", &SolverCallback::getMatrixType) .def("assembleMatrix", &SolverCallback::assembleMatrix) .def("assembleLumpedMatrix", &SolverCallback::assembleLumpedMatrix) .def("assembleResidual", [](SolverCallback & self) { self.assembleResidual(); }) .def("predictor", &SolverCallback::predictor) .def("corrector", &SolverCallback::corrector) .def("beforeSolveStep", &SolverCallback::beforeSolveStep) .def("afterSolveStep", &SolverCallback::afterSolveStep) .def_property_readonly("dof_manager", &SolverCallback::getSCDOFManager, py::return_value_policy::reference); py::class_(mod, "InterceptSolverCallback") .def(py::init_alias()); } } // namespace akantu diff --git a/src/model/phase_field/phase_field_model.cc b/src/model/phase_field/phase_field_model.cc index 900171ccd..e1cf9c449 100644 --- a/src/model/phase_field/phase_field_model.cc +++ b/src/model/phase_field/phase_field_model.cc @@ -1,767 +1,767 @@ /** * @file phase_field_model.cc * * @author Mohit Pundir * * @date creation: Tue Sep 04 2018 * @date last modification: Wed Jun 23 2021 * * @brief Implementation of PhaseFieldModel class * * * @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 "phase_field_model.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" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ PhaseFieldModel::PhaseFieldModel(Mesh & mesh, UInt dim, const ID & id, std::shared_ptr dof_manager, ModelType model_type) : Model(mesh, model_type, std::move(dof_manager), dim, id), phasefield_index("phasefield index", id), phasefield_local_numbering("phasefield local numbering", id) { AKANTU_DEBUG_IN(); 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->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::_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), _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(); 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(); } } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::assignPhaseFieldToElements( const ElementTypeMapArray * filter) { for_each_element( mesh, [&](auto && element) { UInt phase_index = (*phasefield_selector)(element); AKANTU_DEBUG_ASSERT( phase_index < 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::_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::_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, UInt index) { AKANTU_DEBUG_IN(); UInt phase_index = this->phasefield_index(type, _not_ghost)(index); UInt phase_loc_num = this->phasefield_local_numbering(type, _not_ghost)(index); Real energy = this->phasefields[phase_index]->getEnergy(type, phase_loc_num); 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; + // 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(); 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); } // 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); } /* -------------------------------------------------------------------------- */ UInt PhaseFieldModel::getNbData(const Array & elements, const SynchronizationTag & tag) const { UInt size = 0; UInt nb_nodes_per_element = 0; for (const Element & el : elements) { nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); } switch (tag) { case SynchronizationTag::_phasefield_id: { size += elements.size() * sizeof(UInt); break; } case SynchronizationTag::_for_dump: { // damage size += nb_nodes_per_element * sizeof(Real); break; } case SynchronizationTag::_pfm_damage: { size += nb_nodes_per_element * sizeof(Real); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } return size; } /* -------------------------------------------------------------------------- */ 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) { AKANTU_DEBUG_IN(); switch (tag) { case SynchronizationTag::_phasefield_id: { for (auto && element : elements) { UInt recv_phase_index; buffer >> recv_phase_index; UInt & phase_index = phasefield_index(element); if (phase_index != UInt(-1)) { continue; } // add ghosts element to the correct phasefield phase_index = recv_phase_index; UInt 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(); } /* -------------------------------------------------------------------------- */ UInt PhaseFieldModel::getNbData(const Array & dofs, const SynchronizationTag & tag) const { UInt size = 0; switch (tag) { case SynchronizationTag::_for_dump: { size += sizeof(Real); break; } case SynchronizationTag::_pfm_damage: { size += sizeof(Real); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } return size; } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const { switch (tag) { case SynchronizationTag::_for_dump: { packDOFDataHelper(*damage, buffer, dofs); break; } case SynchronizationTag::_pfm_damage: { packDOFDataHelper(*damage, buffer, dofs); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) { switch (tag) { case SynchronizationTag::_for_dump: { unpackDOFDataHelper(*damage, buffer, dofs); break; } case SynchronizationTag::_pfm_damage: { unpackDOFDataHelper(*damage, buffer, dofs); 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*/, UInt /*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/phasefield.cc b/src/model/phase_field/phasefield.cc index b7b7bae49..36be6c7ff 100644 --- a/src/model/phase_field/phasefield.cc +++ b/src/model/phase_field/phasefield.cc @@ -1,420 +1,420 @@ /** * @file phasefield.cc * * @author Mohit Pundir * * @date creation: Fri Jun 19 2020 * @date last modification: Fri May 14 2021 * * @brief Implementation of the common part of the phasefield class * * * @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 "phasefield.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), spatial_dimension(this->model.getSpatialDimension()), element_filter("element_filter", id), damage_on_qpoints("damage", *this), phi("phi", *this), strain("strain", *this), gradd("grad_d", *this), driving_force("driving_force", *this), driving_energy("driving_energy", *this), dissipated_energy("dissipated_energy", *this), damage_energy("damage_energy", *this), damage_energy_density("damage_energy_density", *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, UInt dim, const Mesh & mesh, FEEngine & fe_engine, const ID & id) : Parsable(ParserType::_phasefield, id), id(id), fem(fe_engine), model(model), spatial_dimension(this->model.getSpatialDimension()), element_filter("element_filter", id), damage_on_qpoints("damage", *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), gradd("grad_d", *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), dissipated_energy("dissipated_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) { AKANTU_DEBUG_IN(); /// for each connectivity types allocate the element filer array of the /// material element_filter.initialize(mesh, _spatial_dimension = spatial_dimension, _element_kind = _ek_regular); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ PhaseField::~PhaseField() = default; /* -------------------------------------------------------------------------- */ void PhaseField::initialize() { registerParam("name", name, std::string(), _pat_parsable | _pat_readable); registerParam("l0", l0, Real(0.), _pat_parsable | _pat_readable, "length scale parameter"); registerParam("gc", g_c, _pat_parsable | _pat_readable, "critical local fracture energy density"); registerParam("E", E, _pat_parsable | _pat_readable, "Young's modulus"); registerParam("nu", nu, _pat_parsable | _pat_readable, "Poisson ratio"); damage_on_qpoints.initialize(1); phi.initialize(1); driving_force.initialize(1); driving_energy.initialize(spatial_dimension); gradd.initialize(spatial_dimension); 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_uint.begin(); it != internal_vectors_uint.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(); UInt 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(); const Array & damage = model.getDamage(); 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); // damage_energy_density_on_qpoints = gc/l0 + phi = scalar //auto & damage_energy_density_vect = damage_energy_density(type, ghost_type); auto & driving_force_vect = driving_force(type, ghost_type); Array nt_driving_force(nb_quadrature_points, nb_nodes_per_element); fem.computeNtb(driving_force_vect, nt_driving_force, type, ghost_type, elem_filter); Array int_nt_driving_force(nb_element * nb_quadrature_points, nb_nodes_per_element); fem.integrate(nt_driving_force, int_nt_driving_force, nb_nodes_per_element, type, ghost_type, elem_filter); // damage_energy_on_qpoints = gc*l0 = scalar auto & driving_energy_vect = driving_energy(type, ghost_type); Array bt_driving_energy(nb_element * nb_quadrature_points, nb_nodes_per_element); fem.computeBtD(driving_energy_vect, bt_driving_energy, type, ghost_type, elem_filter); Array int_bt_driving_energy(nb_element, nb_nodes_per_element); fem.integrate(bt_driving_energy, int_bt_driving_energy, nb_nodes_per_element, type, ghost_type, elem_filter); model.getDOFManager().assembleElementalArrayLocalArray( int_nt_driving_force, internal_force, type, ghost_type, -1, elem_filter); model.getDOFManager().assembleElementalArrayLocalArray( int_bt_driving_energy, internal_force, type, ghost_type, -1, elem_filter); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseField::assembleStiffnessMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix"); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { auto & elem_filter = element_filter(type, ghost_type); if (elem_filter.empty()) { AKANTU_DEBUG_OUT(); return; } auto nb_element = elem_filter.size(); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); auto nt_b_n = std::make_unique>( nb_element * nb_quadrature_points, nb_nodes_per_element * nb_nodes_per_element, "N^t*b*N"); auto bt_d_b = std::make_unique>( nb_element * nb_quadrature_points, nb_nodes_per_element * nb_nodes_per_element, "B^t*D*B"); // damage_energy_density_on_qpoints = gc/l0 + phi = scalar auto & damage_energy_density_vect = damage_energy_density(type, ghost_type); // damage_energy_on_qpoints = gc*l0 = scalar auto & damage_energy_vect = damage_energy(type, ghost_type); fem.computeBtDB(damage_energy_vect, *bt_d_b, 2, type, ghost_type, elem_filter); fem.computeNtbN(damage_energy_density_vect, *nt_b_n, type, ghost_type, elem_filter); /// compute @f$ K_{\grad d} = \int_e \mathbf{N}^t * \mathbf{w} * /// \mathbf{N}@f$ auto K_n = std::make_unique>( nb_element, nb_nodes_per_element * nb_nodes_per_element, "K_n"); fem.integrate(*nt_b_n, *K_n, nb_nodes_per_element * nb_nodes_per_element, type, ghost_type, elem_filter); model.getDOFManager().assembleElementalMatricesToMatrix( "K", "damage", *K_n, type, _not_ghost, _symmetric, elem_filter); /// compute @f$ K_{\grad d} = \int_e \mathbf{B}^t * \mathbf{W} * /// \mathbf{B}@f$ auto K_b = std::make_unique>( nb_element, nb_nodes_per_element * nb_nodes_per_element, "K_b"); fem.integrate(*bt_d_b, *K_b, nb_nodes_per_element * nb_nodes_per_element, type, ghost_type, elem_filter); model.getDOFManager().assembleElementalMatricesToMatrix( "K", "damage", *K_b, type, _not_ghost, _symmetric, elem_filter); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseField::computeDissipatedEnergyByElements() { AKANTU_DEBUG_IN(); const Array & damage = model.getDamage(); for (auto type : element_filter.elementTypes(spatial_dimension, _not_ghost)) { Array & elem_filter = element_filter(type, _not_ghost); if (elem_filter.empty()) { continue; } Array & damage_interpolated = damage_on_qpoints(type, _not_ghost); // compute the damage on quadrature points fem.interpolateOnIntegrationPoints(damage, damage_interpolated, 1, type, _not_ghost); Array & gradd_vect = gradd(type, _not_ghost); /// compute @f$\nabla u@f$ fem.gradientOnIntegrationPoints(damage, gradd_vect, 1, type, _not_ghost, elem_filter); computeDissipatedEnergy(type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseField::computeDissipatedEnergy(ElementType /*unused*/) { AKANTU_DEBUG_IN(); AKANTU_TO_IMPLEMENT(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real PhaseField::getEnergy(){ AKANTU_DEBUG_IN(); Real edis = 0.; computeDissipatedEnergyByElements(); /// integrate the dissipated energy for each type of elements for (auto type : element_filter.elementTypes(spatial_dimension, _not_ghost)) { edis += fem.integrate(dissipated_energy(type, _not_ghost), type, _not_ghost, element_filter(type, _not_ghost)); } AKANTU_DEBUG_OUT(); return edis; } /* -------------------------------------------------------------------------- */ Real PhaseField::getEnergy(ElementType type, UInt 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, type, element_filter(type)(index)); return edis; } /* -------------------------------------------------------------------------- */ void PhaseField::beforeSolveStep() { this->savePreviousState(); - //this->computeAllDrivingForces(_not_ghost); + 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/phasefields/phasefield_exponential.hh b/src/model/phase_field/phasefields/phasefield_exponential.hh index c8774c35f..ea291c7cf 100644 --- a/src/model/phase_field/phasefields/phasefield_exponential.hh +++ b/src/model/phase_field/phasefields/phasefield_exponential.hh @@ -1,168 +1,168 @@ /** * @file phasefield_exponential.hh * * @author Mohit Pundir * * @date creation: Fri Jun 19 2020 * @date last modification: Wed Jun 23 2021 * * @brief Phasefield law for approximating discrete crack as an exponential * * * @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 "phasefield.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_PHASEFIELD_EXPONENTIAL_HH__ #define __AKANTU_PHASEFIELD_EXPONENTIAL_HH__ namespace akantu { class PhaseFieldExponential : public PhaseField { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: PhaseFieldExponential(PhaseFieldModel & model, const ID & id = ""); ~PhaseFieldExponential() override = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public : /// compute the dissiapted energy void computeDissipatedEnergy(ElementType el_type) override; void computeDissipatedEnergyByElement(ElementType type, UInt index, Vector & epot_on_quad_points) override; protected: void computePhiOnQuad(const Matrix & /*strain_quad*/, Real & /*phi_quad*/, Real & /*phi_hist_quad*/); void computeDrivingForce(const ElementType & /*el_type*/, GhostType /*ghost_type*/) override; inline void computeDrivingForceOnQuad(const Real & /*phi_quad*/, Real & /*driving_force_quad*/); inline void computeDamageEnergyDensityOnQuad(const Real & /*phi_quad*/, Real & /*dam_energy_quad*/); inline void computeDissipatedEnergyOnQuad(const Real & /*dam_quad*/, const Vector & /*grad_d_quad */, Real & /*energy*/); public: void updateInternalParameters() override; }; /* -------------------------------------------------------------------------- */ inline void PhaseFieldExponential::computeDissipatedEnergyOnQuad( const Real & dam, const Vector & grad_d, Real & edis) { for (auto i : arange(spatial_dimension)) { edis = this->l0*grad_d[i]*grad_d[i]; } edis += this->g_c*dam*dam/(4*this->l0); } /* -------------------------------------------------------------------------- */ inline void PhaseFieldExponential::computeDrivingForceOnQuad(const Real & phi_quad, Real & driving_force_quad) { driving_force_quad = 2.0 * phi_quad; } /* -------------------------------------------------------------------------- */ inline void PhaseFieldExponential::computeDamageEnergyDensityOnQuad( const Real & phi_quad, Real & dam_energy_quad) { dam_energy_quad = 2.0 * phi_quad + this->g_c / this->l0; } /* -------------------------------------------------------------------------- */ inline void PhaseFieldExponential::computePhiOnQuad(const Matrix & strain_quad, Real & phi_quad, Real & phi_hist_quad) { Matrix strain_plus(spatial_dimension, spatial_dimension); Matrix strain_minus(spatial_dimension, spatial_dimension); Matrix strain_dir(spatial_dimension, spatial_dimension); Matrix strain_diag_plus(spatial_dimension, spatial_dimension); Matrix strain_diag_minus(spatial_dimension, spatial_dimension); Vector strain_values(spatial_dimension); Real trace_plus, trace_minus; strain_plus.zero(); strain_minus.zero(); strain_dir.zero(); strain_values.zero(); strain_diag_plus.zero(); strain_diag_minus.zero(); strain_quad.eig(strain_values, strain_dir); for (UInt i = 0; i < spatial_dimension; i++) { strain_diag_plus(i, i) = std::max(Real(0.), strain_values(i)); strain_diag_minus(i, i) = std::min(Real(0.), strain_values(i)); } Matrix mat_tmp(spatial_dimension, spatial_dimension); Matrix sigma_plus(spatial_dimension, spatial_dimension); Matrix sigma_minus(spatial_dimension, spatial_dimension); mat_tmp.mul(strain_diag_plus, strain_dir); strain_plus.mul(strain_dir, mat_tmp); mat_tmp.mul(strain_diag_minus, strain_dir); strain_minus.mul(strain_dir, mat_tmp); trace_plus = std::max(Real(0.), strain_quad.trace()); trace_minus = std::min(Real(0.), strain_quad.trace()); for (UInt i = 0; i < spatial_dimension; i++) { for (UInt j = 0; j < spatial_dimension; j++) { sigma_plus(i, j) = static_cast(i == j) * lambda * trace_plus + 2 * mu * strain_plus(i, j); sigma_minus(i, j) = static_cast(i == j) * lambda * trace_minus + 2 * mu * strain_minus(i, j); } } - phi_quad = 0.5 * sigma_plus.doubleDot(strain_quad); + phi_quad = 0.5 * sigma_plus.doubleDot(strain_plus); if (phi_quad < phi_hist_quad) { phi_quad = phi_hist_quad; } } } // namespace akantu #endif