diff --git a/python/py_dof_manager.cc b/python/py_dof_manager.cc index 6af11fe89..1adddf59d 100644 --- a/python/py_dof_manager.cc +++ b/python/py_dof_manager.cc @@ -1,295 +1,296 @@ /* * 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( "getSolution", [](DOFManager & self) -> decltype(auto) { return self.getSolution(); }, 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("assembleMatMulVectToGlobalArray", &DOFManager::assembleMatMulVectToGlobalArray, py::arg("dof_id"), py::arg("A_id"), py::arg("x"), py::arg("array"), 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) .def("globalToLocalEquationNumber", &DOFManager::globalToLocalEquationNumber) .def("localToGlobalEquationNumber", &DOFManager::localToGlobalEquationNumber) .def("getLocalEquationsNumbers", &DOFManager::getLocalEquationsNumbers) .def("getGlobalBlockedDOFs", &DOFManager::getGlobalBlockedDOFs) .def("updateGlobalBlockedDofs", &DOFManager::updateGlobalBlockedDofs) .def("getTimeStepSolver", &DOFManager::getTimeStepSolver, py::return_value_policy::reference) .def( "assembleElementalArrayLocalArray", [](DOFManager & self, const Array & elementary_vect, Array & array_assembeled, ElementType type, GhostType ghost_type, Real scale_factor) { self.assembleElementalArrayLocalArray(elementary_vect, array_assembeled, type, ghost_type, scale_factor); }, py::arg("elementary_vect"), py::arg("array_assembeled"), py::arg("type"), py::arg("ghost_type") = _not_ghost, py::arg("scale_factor") = 1.) .def( "assembleToGlobalArray", [](DOFManager & self, const ID & dof_id, const Array & array_to_assemble, SolverVector & global_array, Real scale_factor) { self.assembleToGlobalArray(dof_id, array_to_assemble, global_array, scale_factor); }, py::arg("dof_id"), py::arg("array_to_assemble"), py::arg("global_array"), py::arg("scale_factor") = 1.) .def("getNewGlobalVector", &DOFManager::getNewGlobalVector, py::return_value_policy::reference); 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); }) .def("getNbIterations", [](NonLinearSolver & self) -> Int { return self.get("nb_iterations"); }) .def("getConvergenceStatus", [](NonLinearSolver & self) -> bool { return self.get("converged"); }) .def("getError", [](NonLinearSolver & self) -> Real { return self.get("error"); }) .def("getMaxIterations", [](NonLinearSolver & self) -> Int { return self.get("max_iterations"); }) .def("getThreshold", [](NonLinearSolver & self) -> Real { - return self.get("threshold"); + return self.get("threshold_normalized"); }); py::class_(mod, "TimeStepSolver") .def( "assembleResidual", [](TimeStepSolver & self, const ID & residual_part) { self.assembleResidual(residual_part); }, py::arg("residual_part")) .def("assembleResidual", [](TimeStepSolver & self) { self.assembleResidual(); }) + .def("assembleMatrix", &TimeStepSolver::assembleMatrix) .def( "assembleResidual", [](TimeStepSolver & self, SolverCallback & solver_callback, const ID & residual_part) { self.assembleResidual(solver_callback, residual_part); }, py::arg("solver_callback"), py::arg("residual_part")) .def( "assembleResidual", [](TimeStepSolver & self, SolverCallback & solver_callback) { self.assembleResidual(solver_callback); }, py::arg("solver_callback")) .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/python/py_heat_transfer_interface_model.cc b/python/py_heat_transfer_interface_model.cc index c97c07340..a96b0d94c 100644 --- a/python/py_heat_transfer_interface_model.cc +++ b/python/py_heat_transfer_interface_model.cc @@ -1,102 +1,104 @@ /** * @file py_heat_transfer_interface_model.cc * * @author Emil Gallyamov * * @date creation: Tue May 16 2023 * @date last modification: Tue May 16 2023 * * @brief pybind11 interface to HeatTransferInterfaceModel * * * @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 "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ // #include #include // #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ #define def_deprecated(func_name, mesg) \ def(func_name, [](py::args, py::kwargs) { AKANTU_ERROR(mesg); }) #define def_function_nocopy(func_name) \ def( \ #func_name, \ [](HeatTransferInterfaceModel & self) -> decltype(auto) { \ return self.func_name(); \ }, \ py::return_value_policy::reference) #define def_function(func_name) \ def(#func_name, [](HeatTransferInterfaceModel & self) -> decltype(auto) { \ return self.func_name(); \ }) /* -------------------------------------------------------------------------- */ void register_heat_transfer_interface_model(py::module & mod) { py::class_( mod, "HeatTransferInterfaceModel") .def(py::init(), py::arg("mesh"), py::arg("spatial_dimension") = _all_dimensions, py::arg("id") = "heat_transfer_model") .def( "initFull", [](HeatTransferInterfaceModel & self, const HeatTransferModelOptions & options) { self.initFull(options); }, py::arg("_analysis_method") = HeatTransferModelOptions()) .def( "initFull", [](HeatTransferInterfaceModel & self, const AnalysisMethod & _analysis_method) { self.initFull(HeatTransferModelOptions(_analysis_method)); }, py::arg("_analysis_method")) .def("setTimeStep", &HeatTransferInterfaceModel::setTimeStep, py::arg("time_step"), py::arg("solver_id") = "") - .def("getTransversalConductivityOnQpoints", - &HeatTransferInterfaceModel::getTransversalConductivityOnQpoints, - py::arg("el_type"), py::arg("ghost_type") = _not_ghost, - py::return_value_policy::reference) + .def("computeTempOnQpoints", + &HeatTransferInterfaceModel::computeTempOnQpoints, + py::arg("ghost_type") = _not_ghost) .def("getLongitudinalConductivityOnQpoints", &HeatTransferInterfaceModel::getLongitudinalConductivityOnQpoints, py::arg("el_type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def( "getOpening", [](HeatTransferInterfaceModel & self, ElementType & el_type, GhostType & ghost_type) { self.getOpening(el_type, ghost_type); }, py::arg("el_type"), py::arg("ghost_type") = _not_ghost, - py::return_value_policy::reference); + py::return_value_policy::reference) + .def("updateNormalOpeningAtQuadraturePoints", + &HeatTransferInterfaceModel::updateNormalOpeningAtQuadraturePoints, + py::arg("positions"), py::arg("ghost_type") = _not_ghost); } } // namespace akantu diff --git a/python/py_heat_transfer_model.cc b/python/py_heat_transfer_model.cc index 568397bcb..00cd5c9c2 100644 --- a/python/py_heat_transfer_model.cc +++ b/python/py_heat_transfer_model.cc @@ -1,142 +1,143 @@ /** * @file py_heat_transfer_model.cc * * @author Nicolas Richart * * @date creation: Sun Jun 16 2019 * @date last modification: Sun Jun 16 2019 * * @brief pybind11 interface to HeatTransferModel * * * @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 "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ // #include #include // #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ #define def_deprecated(func_name, mesg) \ def(func_name, [](py::args, py::kwargs) { AKANTU_ERROR(mesg); }) #define def_function_nocopy(func_name) \ def( \ #func_name, \ [](HeatTransferModel & self) -> decltype(auto) { \ return self.func_name(); \ }, \ py::return_value_policy::reference) #define def_function(func_name) \ def(#func_name, [](HeatTransferModel & self) -> decltype(auto) { \ return self.func_name(); \ }) /* -------------------------------------------------------------------------- */ void register_heat_transfer_model(py::module & mod) { py::class_(mod, "HeatTransferModelOptions") .def(py::init(), py::arg("analysis_method") = _explicit_lumped_mass); py::class_(mod, "HeatTransferModel", py::multiple_inheritance()) .def(py::init(), py::arg("mesh"), py::arg("spatial_dimension") = _all_dimensions, py::arg("id") = "heat_transfer_model") .def( "initFull", [](HeatTransferModel & self, const HeatTransferModelOptions & options) { self.initFull(options); }, py::arg("_analysis_method") = HeatTransferModelOptions()) .def( "initFull", [](HeatTransferModel & self, const AnalysisMethod & _analysis_method) { self.initFull(HeatTransferModelOptions(_analysis_method)); }, py::arg("_analysis_method")) .def("setTimeStep", &HeatTransferModel::setTimeStep, py::arg("time_step"), py::arg("solver_id") = "") .def("getTimeStep", &HeatTransferModel::getTimeStep, py::arg("solver_id") = "") .def_function(getStableTimeStep) .def_function(getCapacity) .def_function(getDensity) .def_function(assembleInternalHeatRate) + .def_function(assembleCapacity) .def_function_nocopy(getTemperature) .def_function_nocopy(getTemperatureRate) .def_function_nocopy(getBlockedDOFs) .def_function_nocopy(getExternalHeatRate) .def_function_nocopy(getInternalHeatRate) .def_function_nocopy(getMesh) .def_function(needToReassembleCapacity) .def_function(getConductivityMatrixRelease) .def("getConductivityRelease", &HeatTransferModel::getConductivityRelease, py::arg("ghost_type") = _not_ghost) .def("getTemperatureGradient", &HeatTransferModel::getTemperatureGradient, py::arg("el_type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def("getKgradT", &HeatTransferModel::getKgradT, py::arg("el_type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def("getTemperatureRateOnQpoints", &HeatTransferModel::getTemperatureRateOnQpoints, py::arg("el_type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def("applyBC", [](HeatTransferModel & self, BC::Dirichlet::DirichletFunctor & func, const std::string & element_group) { self.applyBC(func, element_group); }) .def("synchronizeField", &HeatTransferModel::synchronizeField, py::arg("synchronization_tag")) .def( "assignPropertyToPhysicalGroup", [](HeatTransferModel & self, const std::string & property_name, const std::string & group_name, Real value) { self.assignPropertyToPhysicalGroup(property_name, group_name, value); }, py::arg("property_name"), py::arg("group_name"), py::arg("value")) .def( "assignPropertyToPhysicalGroup", [](HeatTransferModel & self, const std::string & property_name, const std::string & group_name, Matrix cond_matrix) { self.assignPropertyToPhysicalGroup(property_name, group_name, cond_matrix); }, py::arg("property_name"), py::arg("group_name"), py::arg("cond_matrix")); } } // namespace akantu diff --git a/python/py_material.cc b/python/py_material.cc index a496d2f8a..c220cbc3f 100644 --- a/python/py_material.cc +++ b/python/py_material.cc @@ -1,478 +1,521 @@ /** * @file py_material.cc * * @author Guillaume Anciaux * @author Mohit Pundir * @author Nicolas Richart * * @date creation: Thu Jun 20 2019 * @date last modification: Fri Apr 09 2021 * * @brief pybind11 interface to Material * * * @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 "py_aka_array.hh" #include "py_akantu_pybind11_compatibility.hh" /* -------------------------------------------------------------------------- */ #include #include #if defined(AKANTU_COHESIVE_ELEMENT) #include #include #include #include #include #endif #include /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { namespace { template class PyMaterial : public _Material { public: /* Inherit the constructors */ using _Material::_Material; ~PyMaterial() override = default; void initMaterial() override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, _Material, initMaterial, ); }; void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE_PURE(void, _Material, computeStress, el_type, ghost_type); } void computeTangentModuli(ElementType el_type, Array & tangent_matrix, GhostType ghost_type = _not_ghost) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, _Material, computeTangentModuli, el_type, tangent_matrix, ghost_type); } - void computePotentialEnergy(ElementType el_type) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, _Material, computePotentialEnergy, el_type); } Real getPushWaveSpeed(const Element & element) const override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(Real, _Material, getPushWaveSpeed, element); } Real getShearWaveSpeed(const Element & element) const override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(Real, _Material, getShearWaveSpeed, element); } template void registerInternal(const std::string & name, UInt nb_component) { auto && internal = std::make_shared>(name, *this); AKANTU_DEBUG_INFO("alloc internal " << name << " " << &this->internals[name]); internal->initialize(nb_component); this->internals[name] = internal; } template void setDefaultValueToInternal(const ID & int_id, const T value) { auto & internal_field = this->template getInternal(int_id); internal_field.setDefaultValue(value); } protected: std::map> internals; }; /* ------------------------------------------------------------------------ */ template void register_material_classes(py::module & mod, const std::string & name) { py::class_<_Material, Material, Parsable, PyMaterial<_Material>>( mod, name.c_str(), py::multiple_inheritance()) .def(py::init()) .def("registerInternalReal", [](Material & self, const std::string & name, Int nb_component) { return dynamic_cast &>(self) .template registerInternal(name, nb_component); }) .def("registerInternalUInt", [](Material & self, const std::string & name, Int nb_component) { return dynamic_cast &>(self) .template registerInternal(name, nb_component); }) .def("setDefaultValueToInternalReal", [](Material & self, const ID & int_id, const Real value) { return dynamic_cast &>(self) .template setDefaultValueToInternal(int_id, value); }) .def("setDefaultValueToInternalUInt", [](Material & self, const ID & int_id, const UInt value) { return dynamic_cast &>(self) .template setDefaultValueToInternal(int_id, value); }); } #if defined(AKANTU_COHESIVE_ELEMENT) // trampoline for the cohesive materials template class PyMaterialCohesive : public PyMaterial<_Material> { using Parent = PyMaterial<_Material>; public: using Parent::Parent; void checkInsertion(bool check_only) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, _Material, checkInsertion, check_only); } void computeTraction(const Array & normal, ElementType el_type, GhostType ghost_type = _not_ghost) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE_PURE(void, _Material, computeTraction, normal, el_type, ghost_type); } void computeTangentTraction(ElementType el_type, Array & tangent_matrix, const Array & normal, GhostType ghost_type = _not_ghost) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, _Material, computeTangentTraction, el_type, tangent_matrix, normal, ghost_type); } void computeStress(ElementType /*el_type*/, GhostType /*ghost_type*/ = _not_ghost) final {} void computeTangentModuli(ElementType /*el_type*/, Array & /*tangent_matrix*/, GhostType /*ghost_type*/ = _not_ghost) final {} template void registerInternal(const std::string & name, Int nb_component) { auto && internal = std::make_shared>(name, *this); AKANTU_DEBUG_INFO("alloc internal " << name << " " << &this->internals[name]); internal->initialize(nb_component); this->internals[name] = internal; } }; // trampoline for the cohesive material inheritance where computeTraction is // not pure virtual template class PyMaterialCohesiveDaughters : public PyMaterialCohesive<_Material> { using Parent = PyMaterialCohesive<_Material>; public: using Parent::Parent; void computeTraction(const Array & normal, ElementType el_type, GhostType ghost_type = _not_ghost) override { // NOLINTNEXTLINE PYBIND11_OVERRIDE(void, _Material, computeTraction, normal, el_type, ghost_type); } + + void computeTangentTraction(ElementType el_type, + Array & tangent_matrix, + const Array & normal, + GhostType ghost_type = _not_ghost) override { + // NOLINTNEXTLINE + PYBIND11_OVERRIDE(void, _Material, computeTangentTraction, el_type, + tangent_matrix, normal, ghost_type); + } }; template void register_material_cohesive_classes(py::module & mod, const std::string & name) { py::class_<_Material, MaterialCohesive, PyMaterialCohesiveDaughters<_Material>>( mod, name.c_str(), py::multiple_inheritance()) .def(py::init()) .def("registerInternalReal", [](_Material & self, const std::string & name, UInt nb_component) { return dynamic_cast &>(self) .template registerInternal(name, nb_component); }) .def("registerInternalUInt", [](_Material & self, const std::string & name, UInt nb_component) { return dynamic_cast &>(self) .template registerInternal(name, nb_component); }) .def("setDefaultValueToInternalReal", [](_Material & self, const ID & int_id, const Real value) { return dynamic_cast &>(self) .template setDefaultValueToInternal(int_id, value); }) .def("setDefaultValueToInternalUInt", [](_Material & self, const ID & int_id, const UInt value) { return dynamic_cast &>(self) .template setDefaultValueToInternal(int_id, value); }) - .def("computeTraction", [](_Material & self, const Array & normal, - ElementType el_type, GhostType ghost_type) { + .def("computeTraction", + [](_Material & self, const Array & normal, + ElementType el_type, GhostType ghost_type) { + return dynamic_cast &>( + self) + .computeTraction(normal, el_type, ghost_type); + }) + .def("computeTangentTraction", [](_Material & self, ElementType el_type, + Array & tangent_matrix, + const Array & normal, + GhostType ghost_type) { return dynamic_cast &>(self) - .computeTraction(normal, el_type, ghost_type); + .computeTangentTraction(el_type, tangent_matrix, normal, + ghost_type); }); ; } #endif /* ------------------------------------------------------------------------ */ template void register_internal_field(py::module & mod, const std::string & name) { py::class_, ElementTypeMapArray, std::shared_ptr>>( mod, ("InternalField" + name).c_str()); } } // namespace /* -------------------------------------------------------------------------- */ void register_material(py::module & mod) { py::class_(mod, "MaterialFactory") .def_static( "getInstance", []() -> MaterialFactory & { return Material::getFactory(); }, py::return_value_policy::reference) .def("registerAllocator", [](MaterialFactory & self, const std::string id, py::function func) { self.registerAllocator( id, [func, id](UInt dim, const ID & /*unused*/, SolidMechanicsModel & model, const ID & option) -> std::unique_ptr { py::object obj = func(dim, id, model, option); auto & ptr = py::cast(obj); obj.release(); return std::unique_ptr(&ptr); }); }) .def("getPossibleAllocators", &MaterialFactory::getPossibleAllocators); register_internal_field(mod, "Real"); register_internal_field(mod, "UInt"); py::class_>( mod, "Material", py::multiple_inheritance()) .def(py::init()) .def( "getGradU", [](Material & self, ElementType el_type, GhostType ghost_type = _not_ghost) -> decltype(auto) { return self.getGradU(el_type, ghost_type); }, py::arg("el_type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def( "getStress", [](Material & self, ElementType el_type, GhostType ghost_type = _not_ghost) -> decltype(auto) { return self.getStress(el_type, ghost_type); }, py::arg("el_type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def( "getPotentialEnergy", [](Material & self, ElementType el_type) -> decltype(auto) { return self.getPotentialEnergy(el_type); }, py::arg("el_type"), py::return_value_policy::reference) .def( "getPotentialEnergy", [](Material & self, ElementType el_type, UInt index) -> Real { return self.getPotentialEnergy(el_type, index); }, py::arg("el_type"), py::arg("index")) .def("getPotentialEnergy", [](Material & self) -> Real { return self.getPotentialEnergy(); }) .def("initMaterial", &Material::initMaterial) .def("getModel", &Material::getModel) .def("registerInternalReal", [](Material & self, const std::string & name, UInt nb_component) { return dynamic_cast &>(self) .registerInternal(name, nb_component); }) .def("registerInternalUInt", [](Material & self, const std::string & name, UInt nb_component) { return dynamic_cast &>(self) .registerInternal(name, nb_component); }) .def("setDefaultValueToInternalReal", [](Material & self, const ID & int_id, const Real value) { return dynamic_cast &>(self) .setDefaultValueToInternal(int_id, value); }) .def("setDefaultValueToInternalUInt", [](Material & self, const ID & int_id, const UInt value) { return dynamic_cast &>(self) .setDefaultValueToInternal(int_id, value); }) .def( "getInternalReal", [](Material & self, const ID & id) -> decltype(auto) { return self.getInternal(id); }, py::arg("id"), py::return_value_policy::reference) + .def( + "getInternalRealPrevious", + [](Material & self, const ID & id) -> decltype(auto) { + return self.getInternal(id).previous(); + }, + py::arg("id"), py::return_value_policy::reference) .def( "getInternalUInt", [](Material & self, const ID & id) -> decltype(auto) { return self.getInternal(id); }, py::arg("id"), py::return_value_policy::reference) .def( "getElementFilter", [](Material & self) -> decltype(auto) { return self.getElementFilter(); }, py::return_value_policy::reference) /* * These functions override the `Parsable` interface. * This ensure that the `updateInternalParameters()` function is called. */ .def( "setReal", [](Material & self, const ID & name, const Real value) -> void { self.setParam(name, value); return; }, py::arg("name"), py::arg("value")) .def( "setBool", [](Material & self, const ID & name, const bool value) -> void { self.setParam(name, value); return; }, py::arg("name"), py::arg("value")) .def( "setString", [](Material & self, const ID & name, const std::string & value) -> void { self.setParam(name, value); return; }, py::arg("name"), py::arg("value")) .def( "setInt", [](Material & self, const ID & name, const int value) -> void { self.setParam(name, value); return; }, py::arg("name"), py::arg("value")) .def("getPushWaveSpeed", &Material::getPushWaveSpeed) .def("getShearWaveSpeed", &Material::getShearWaveSpeed) .def("__repr__", [](Material & self) { std::stringstream sstr; sstr << self; return sstr.str(); }); register_material_classes>(mod, "MaterialElastic2D"); register_material_classes>(mod, "MaterialElastic3D"); #if defined(AKANTU_COHESIVE_ELEMENT) /* ------------------------------------------------------------------------ */ py::class_>(mod, "MaterialCohesive", py::multiple_inheritance()) .def(py::init()) .def("registerInternalReal", [](MaterialCohesive & self, const std::string & name, UInt nb_component) { return dynamic_cast &>(self) .registerInternal(name, nb_component); }) .def("registerInternalUInt", [](MaterialCohesive & self, const std::string & name, UInt nb_component) { return dynamic_cast &>(self) .registerInternal(name, nb_component); }) .def( "getFacetFilter", [](MaterialCohesive & self) -> decltype(auto) { return self.getFacetFilter(); }, py::return_value_policy::reference) .def( "getFacetFilter", [](MaterialCohesive & self, ElementType type, GhostType ghost_type) -> decltype(auto) { return self.getFacetFilter(type, ghost_type); }, py::arg("type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def( "getOpening", [](MaterialCohesive & self, ElementType type, GhostType ghost_type) -> decltype(auto) { return self.getOpening(type, ghost_type); }, py::arg("type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def( "getTraction", [](MaterialCohesive & self, ElementType type, GhostType ghost_type) -> decltype(auto) { return self.getTraction(type, ghost_type); }, py::arg("type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def( "computeTraction", [](MaterialCohesive & self, GhostType ghost_type) { return self.computeTraction(ghost_type); }, py::arg("ghost_type") = _not_ghost) .def( "computeTraction", [](MaterialCohesive & self, const Array & normal, ElementType el_type, GhostType ghost_type) { return dynamic_cast &>(self) .computeTraction(normal, el_type, ghost_type); }, py::arg("normal"), py::arg("el_type"), py::arg("ghost_type") = _not_ghost) + .def( + "computeTangentTraction", + [](MaterialCohesive & self, ElementType el_type, + Array & tangent_matrix, const Array & normal, + GhostType ghost_type) { + return dynamic_cast &>(self) + .computeTangentTraction(el_type, tangent_matrix, normal, + ghost_type); + }, + py::arg("el_type"), py::arg("tangent_matrix"), py::arg("normal"), + py::arg("ghost_type") = _not_ghost) + .def( + "computeNormal", + [](MaterialCohesive & self, const Array & position, + Array & normal, ElementType type, GhostType ghost_type) { + self.computeNormal(position, normal, type, ghost_type); + }, + py::arg("position"), py::arg("normal"), py::arg("type"), + py::arg("ghost_type") = _not_ghost) .def("getNormalsAtQuads", &MaterialCohesive::getNormalsAtQuads); register_material_cohesive_classes>( mod, "MaterialCohesiveLinear2D"); register_material_cohesive_classes>( mod, "MaterialCohesiveLinear3D"); register_material_cohesive_classes>( mod, "MaterialCohesiveBilinear2D"); register_material_cohesive_classes>( mod, "MaterialCohesiveBilinear3D"); register_material_cohesive_classes>( mod, "MaterialCohesiveLinearFriction2D"); register_material_cohesive_classes>( mod, "MaterialCohesiveLinearFriction3D"); #endif } } // namespace akantu diff --git a/python/py_mesh.cc b/python/py_mesh.cc index e379a189c..140d86288 100644 --- a/python/py_mesh.cc +++ b/python/py_mesh.cc @@ -1,290 +1,292 @@ /** * @file py_mesh.cc * * @author Guillaume Anciaux * @author Philip Mueller * @author Mohit Pundir * @author Nicolas Richart * * @date creation: Sun Jun 16 2019 * @date last modification: Mon Mar 15 2021 * * @brief pybind11 interface to Mesh * * * @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 "aka_config.hh" /* -------------------------------------------------------------------------- */ #include "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { namespace { /* ------------------------------------------------------------------------ */ template decltype(auto) register_element_type_map_array(py::module & mod, const std::string & name) { return py::class_, std::shared_ptr>>( mod, ("ElementTypeMapArray" + name).c_str()) .def(py::init(), py::arg("id") = "by_element_type_array", py::arg("parent_id") = "no_parent") .def( "__call__", [](ElementTypeMapArray & self, ElementType type, GhostType ghost_type) -> decltype(auto) { return self(type, ghost_type); }, py::arg("type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference, py::keep_alive<0, 1>()) .def( "elementTypes", [](ElementTypeMapArray & self, UInt _dim, GhostType _ghost_type, ElementKind _kind) -> std::vector { auto types = self.elementTypes(_dim, _ghost_type, _kind); std::vector _types; for (auto && t : types) { _types.push_back(t); } return _types; }, py::arg("dim") = _all_dimensions, py::arg("ghost_type") = _not_ghost, py::arg("kind") = _ek_regular) .def( "initialize", [](ElementTypeMapArray & self, const Mesh & mesh, GhostType ghost_type = _casper, UInt nb_component = 1, UInt spatial_dimension = UInt(-2), ElementKind element_kind = _ek_not_defined, bool with_nb_element = false, bool with_nb_nodes_per_element = false, T default_value = T(), bool do_not_default = false) { self.initialize( mesh, _ghost_type = ghost_type, _nb_component = nb_component, _spatial_dimension = (spatial_dimension == UInt(-2) ? mesh.getSpatialDimension() : spatial_dimension), _element_kind = element_kind, _with_nb_element = with_nb_element, _with_nb_nodes_per_element = with_nb_nodes_per_element, _default_value = default_value, _do_not_default = do_not_default); }, py::arg("mesh"), py::arg("ghost_type") = _casper, py::arg("nb_component") = 1, py::arg("spatial_dimension") = UInt(-2), py::arg("element_kind") = _ek_not_defined, py::arg("with_nb_element") = false, py::arg("with_nb_nodes_per_element") = false, py::arg("default_value") = T(), py::arg("do_not_default") = false); } } // namespace /* -------------------------------------------------------------------------- */ void register_mesh(py::module & mod) { py::class_(mod, "PeriodicSlaves") .def( "__iter__", [](Mesh::PeriodicSlaves & _this) { py::make_iterator(_this.begin(), _this.end()); }, py::keep_alive<0, 1>()); py::class_(mod, "MeshData") .def( "getElementalDataUInt", [](MeshData & _this, const ID & name) -> decltype(auto) { return _this.getElementalData(name); }, py::return_value_policy::reference) .def( "getElementalDataReal", [](MeshData & _this, const ID & name) -> decltype(auto) { return _this.getElementalData(name); }, py::return_value_policy::reference); py::class_(mod, "Mesh", py::multiple_inheritance()) .def(py::init(), py::arg("spatial_dimension"), py::arg("id") = "mesh") .def("read", &Mesh::read, py::arg("filename"), py::arg("mesh_io_type") = _miot_auto, "read the mesh from a file") .def( "getNodes", [](Mesh & self) -> decltype(auto) { return self.getNodes(); }, py::return_value_policy::reference) .def("getNbNodes", &Mesh::getNbNodes) + .def("getLowerBounds", &Mesh::getLowerBounds) + .def("getUpperBounds", &Mesh::getUpperBounds) .def( "getNodesFlags", [](const Mesh & self) -> decltype(auto) { return self.getNodesFlags(); }, py::return_value_policy::reference) .def( "getConnectivity", [](Mesh & self, ElementType type) -> decltype(auto) { return self.getConnectivity(type); }, py::return_value_policy::reference) .def( "getConnectivity", [](Mesh & self, ElementType type, GhostType ghost_type) -> decltype(auto) { return self.getConnectivity(type, ghost_type); }, py::return_value_policy::reference) .def( "getConnectivities", [](Mesh & self) -> decltype(auto) { return self.getConnectivities(); }, py::return_value_policy::reference) .def( "addConnectivityType", [](Mesh & self, ElementType type, GhostType ghost_type) -> void { self.addConnectivityType(type, ghost_type); }, py::arg("type"), py::arg("ghost_type") = _not_ghost) .def("distribute", [](Mesh & self) { self.distribute(); }) .def("isDistributed", [](const Mesh & self) { return self.isDistributed(); }) .def("fillNodesToElements", &Mesh::fillNodesToElements, py::arg("dimension") = _all_dimensions) .def("getAssociatedElements", [](Mesh & self, const UInt & node, py::list list) { Array elements; self.getAssociatedElements(node, elements); for (auto && element : elements) { list.append(element); } }) .def("makePeriodic", [](Mesh & self, const SpatialDirection & direction) { self.makePeriodic(direction); }) .def( "getNbElement", [](Mesh & self, const UInt spatial_dimension, GhostType ghost_type, ElementKind kind) { return self.getNbElement(spatial_dimension, ghost_type, kind); }, py::arg("spatial_dimension") = _all_dimensions, py::arg("ghost_type") = _not_ghost, py::arg("kind") = _ek_not_defined) .def( "getNbElement", [](Mesh & self, ElementType type, GhostType ghost_type) { return self.getNbElement(type, ghost_type); }, py::arg("type"), py::arg("ghost_type") = _not_ghost) .def_static( "getSpatialDimension", [](ElementType & type) { return Mesh::getSpatialDimension(type); }) .def_static("getNaturalSpaceDimension", [](ElementType & type) { return Mesh::getNaturalSpaceDimension(type); }) .def_static("getKind", [](ElementType & type) { return Mesh::getKind(type); }) .def_static( "getNbNodesPerElement", [](ElementType & type) { return Mesh::getNbNodesPerElement(type); }) .def( "getDataReal", [](Mesh & _this, const ID & name, ElementType type, GhostType ghost_type) -> decltype(auto) { return _this.getData(name, type, ghost_type); }, py::arg("name"), py::arg("type"), py::arg("ghost_type") = _not_ghost, py::return_value_policy::reference) .def( "hasDataReal", [](Mesh & _this, const ID & name, ElementType type, GhostType ghost_type) -> bool { return _this.hasData(name, type, ghost_type); }, py::arg("name"), py::arg("type"), py::arg("ghost_type") = _not_ghost) .def("isPeriodic", [](const Mesh & _this) { return _this.isPeriodic(); }) .def("getPeriodicMaster", &Mesh::getPeriodicMaster) .def("getPeriodicSlaves", &Mesh::getPeriodicSlaves) .def("isPeriodicSlave", &Mesh::isPeriodicSlave) .def("isPeriodicMaster", &Mesh::isPeriodicMaster) .def( "getMeshFacets", [](const Mesh & self) -> const Mesh & { return self.getMeshFacets(); }, py::return_value_policy::reference) .def("initMeshFacets", &Mesh::initMeshFacets, py::arg("id") = "mesh_facets", py::return_value_policy::reference) .def( "elementTypes", [](Mesh & self, UInt spatial_dimension, GhostType ghost_type, ElementKind kind) -> std::vector { auto types = self.elementTypes(spatial_dimension, ghost_type, kind); std::vector _types; for (auto && t : types) { _types.push_back(t); } return _types; }, py::arg("spatial_dimension") = _all_dimensions, py::arg("ghost_type") = _not_ghost, py::arg("kind") = _ek_regular); /* ------------------------------------------------------------------------ */ py::class_(mod, "MeshUtils") .def_static("buildFacets", &MeshUtils::buildFacets); py::class_(mod, "MeshAccessor") .def(py::init(), py::arg("mesh")) .def( "resizeConnectivity", [](MeshAccessor & self, UInt new_size, ElementType type, GhostType gt) -> void { self.resizeConnectivity(new_size, type, gt); }, py::arg("new_size"), py::arg("type"), py::arg("ghost_type") = _not_ghost) .def( "resizeNodes", [](MeshAccessor & self, UInt new_size) -> void { self.resizeNodes(new_size); }, py::arg("new_size")) .def("makeReady", &MeshAccessor::makeReady); register_element_type_map_array(mod, "Real"); register_element_type_map_array(mod, "UInt"); register_element_type_map_array(mod, "bool"); // register_element_type_map_array(mod, "String"); } } // namespace akantu diff --git a/python/py_model.cc b/python/py_model.cc index 96b091cb9..deaac9af1 100644 --- a/python/py_model.cc +++ b/python/py_model.cc @@ -1,150 +1,151 @@ /** * @file py_model.cc * * @author Guillaume Anciaux * @author Emil Gallyamov * @author Philip Mueller * @author Mohit Pundir * @author Nicolas Richart * * @date creation: Sun Jun 16 2019 * @date last modification: Sat Mar 13 2021 * * @brief pybind11 interface to Model and parent classes * * * @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 "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include #include /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ void register_model(py::module & mod) { py::class_(mod, "ModelSolver", py::multiple_inheritance()) .def( "getNonLinearSolver", [](ModelSolver & self, const ID & solver_id) -> NonLinearSolver & { return self.getNonLinearSolver(solver_id); }, py::arg("solver_id") = "", py::return_value_policy::reference) .def( "getTimeStepSolver", [](ModelSolver & self, const ID & solver_id) -> TimeStepSolver & { return self.getTimeStepSolver(solver_id); }, py::arg("solver_id") = "", py::return_value_policy::reference) .def( "solveStep", [](ModelSolver & self, const ID & solver_id) { self.solveStep(solver_id); }, py::arg("solver_id") = "") .def( "solveStep", [](ModelSolver & self, SolverCallback & callback, const ID & solver_id) { self.solveStep(callback, solver_id); }, py::arg("callback"), py::arg("solver_id") = ""); py::class_(mod, "Model", py::multiple_inheritance()) .def("getSpatialDimension", &Model::getSpatialDimension) .def("setBaseName", &Model::setBaseName) .def("setDirectory", &Model::setDirectory) .def("getFEEngine", &Model::getFEEngine, py::arg("name") = "", py::return_value_policy::reference) .def("getFEEngineBoundary", &Model::getFEEngine, py::arg("name") = "", py::return_value_policy::reference) .def("addDumpFieldVector", &Model::addDumpFieldVector) .def("addDumpField", &Model::addDumpField) .def("setBaseNameToDumper", &Model::setBaseNameToDumper) + .def("setDirectoryToDumper", &Model::setDirectoryToDumper) .def("addDumpFieldVectorToDumper", &Model::addDumpFieldVectorToDumper) .def("addDumpFieldToDumper", &Model::addDumpFieldToDumper) .def("dump", [](Model & self) { self.dump(); }) .def( "dump", [](Model & self, UInt step) { self.dump(step); }, py::arg("step")) .def( "dump", [](Model & self, Real time, UInt step) { self.dump(time, step); }, py::arg("time"), py::arg("step")) .def( "dump", [](Model & self, const std::string & dumper) { self.dump(dumper); }, py::arg("dumper_name")) .def( "dump", [](Model & self, const std::string & dumper, UInt step) { self.dump(dumper, step); }, py::arg("dumper_name"), py::arg("step")) .def( "dump", [](Model & self, const std::string & dumper, Real time, UInt step) { self.dump(dumper, time, step); }, py::arg("dumper_name"), py::arg("time"), py::arg("step")) .def("initNewSolver", &Model::initNewSolver) .def( "getNewSolver", [](Model & self, const std::string id, const TimeStepSolverType & time, const NonLinearSolverType & type) { self.getNewSolver(id, time, type); }, py::return_value_policy::reference) .def( "setIntegrationScheme", [](Model & self, const std::string id, const std::string primal, const IntegrationSchemeType & scheme_type, IntegrationScheme::SolutionType solution_type) { self.setIntegrationScheme(id, primal, scheme_type, solution_type); }, py::arg("id"), py::arg("primal"), py::arg("scheme_type"), py::arg("solution_type") = IntegrationScheme::SolutionType::_not_defined) // .def("setIntegrationScheme", // [](Model & self, const std::string id, const std::string primal, // std::unique_ptr & scheme, // IntegrationScheme::SolutionType solution_type) { // self.setIntegrationScheme(id, primal, scheme, solution_type); // }) .def("getDOFManager", &Model::getDOFManager, py::return_value_policy::reference) .def("assembleMatrix", &Model::assembleMatrix); } } // namespace akantu diff --git a/python/py_solid_mechanics_model.cc b/python/py_solid_mechanics_model.cc index 2446e2b86..899af20da 100644 --- a/python/py_solid_mechanics_model.cc +++ b/python/py_solid_mechanics_model.cc @@ -1,187 +1,199 @@ /** * @file py_solid_mechanics_model.cc * * @author Guillaume Anciaux * @author Mohit Pundir * @author Nicolas Richart * * @date creation: Sun Jun 16 2019 * @date last modification: Sat Mar 13 2021 * * @brief pybind11 interface to SolidMechanicsModel * * * @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 "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ #define def_deprecated(func_name, mesg) \ def(func_name, [](py::args, py::kwargs) { AKANTU_ERROR(mesg); }) #define def_function_nocopy(func_name) \ def( \ #func_name, \ [](SolidMechanicsModel & self) -> decltype(auto) { \ return self.func_name(); \ }, \ py::return_value_policy::reference) #define def_function(func_name) \ def(#func_name, [](SolidMechanicsModel & self) -> decltype(auto) { \ return self.func_name(); \ }) /* -------------------------------------------------------------------------- */ void register_solid_mechanics_model(py::module & mod) { py::class_(mod, "SolidMechanicsModelOptions") .def(py::init(), py::arg("_analysis_method") = _explicit_lumped_mass); py::class_(mod, "SolidMechanicsModel", py::multiple_inheritance()) .def(py::init, const ModelType>(), py::arg("mesh"), py::arg("spatial_dimension") = _all_dimensions, py::arg("id") = "solid_mechanics_model", py::arg("dof_manager") = nullptr, py::arg("model_type") = ModelType::_solid_mechanics_model) .def( "initFull", [](SolidMechanicsModel & self, const SolidMechanicsModelOptions & options) { self.initFull(options); }, py::arg("option") = SolidMechanicsModelOptions()) .def( "initFull", [](SolidMechanicsModel & self, const AnalysisMethod & analysis_method) { self.initFull(_analysis_method = analysis_method); }, py::arg("_analysis_method")) .def_deprecated("applyDirichletBC", "Deprecated: use applyBC") .def("applyBC", [](SolidMechanicsModel & self, BC::Dirichlet::DirichletFunctor & func, const std::string & element_group) { self.applyBC(func, element_group); }) .def("applyBC", [](SolidMechanicsModel & self, BC::Neumann::NeumannFunctor & func, const std::string & element_group) { self.applyBC(func, element_group); }) .def("setTimeStep", &SolidMechanicsModel::setTimeStep, py::arg("time_step"), py::arg("solver_id") = "") .def( "getEnergy", [](SolidMechanicsModel & self, const std::string & energy_id) { return self.getEnergy(energy_id); }, py::arg("energy_id")) .def( "getEnergy", [](SolidMechanicsModel & self, const std::string & energy_id, const std::string & group_id) { return self.getEnergy(energy_id, group_id); }, py::arg("energy_id"), py::arg("group_id")) - + .def("getDisplacementRelease", + [](SolidMechanicsModel & self) { + return self.getDisplacementRelease(); + }) + .def("setDisplacementRelease", + [](SolidMechanicsModel & self, UInt value) { + return self.setDisplacementRelease(value); + }) + .def("increaseDisplacementRelease", + [](SolidMechanicsModel & self) { + UInt release = self.getDisplacementRelease(); + return self.setDisplacementRelease(release + 1); + }) .def_function(assembleStiffnessMatrix) .def_function(assembleInternalForces) .def_function(assembleMass) .def_function(assembleMassLumped) .def_function(getStableTimeStep) .def_function_nocopy(getExternalForce) .def_function_nocopy(getDisplacement) .def_function_nocopy(getPreviousDisplacement) .def_function_nocopy(getCurrentPosition) .def_function_nocopy(getIncrement) .def_function_nocopy(getMass) .def_function_nocopy(getVelocity) .def_function_nocopy(getAcceleration) .def_function_nocopy(getInternalForce) .def_function_nocopy(getBlockedDOFs) .def_function_nocopy(getMesh) .def( "getMaterial", [](SolidMechanicsModel & self, UInt material_id) -> decltype(auto) { return self.getMaterial(material_id); }, py::arg("material_id"), py::return_value_policy::reference) .def( "getMaterial", [](SolidMechanicsModel & self, const ID & material_name) -> decltype(auto) { return self.getMaterial(material_name); }, py::arg("material_name"), py::return_value_policy::reference) .def( "getMaterial", [](SolidMechanicsModel & self, const Element & element) -> decltype(auto) { return self.getMaterial(element); }, py::arg("element"), py::return_value_policy::reference) .def("getNbMaterials", &SolidMechanicsModel::getNbMaterials) .def("getMaterialIndex", &SolidMechanicsModel::getMaterialIndex) .def("setMaterialSelector", [](SolidMechanicsModel & self, std::shared_ptr material_selector) { std::cout << (*material_selector)(ElementNull) << std::endl; self.setMaterialSelector(material_selector); }) .def("getMaterialSelector", &SolidMechanicsModel::getMaterialSelector) .def( "getMaterialByElement", [](const SolidMechanicsModel & self) -> decltype(auto) { return self.getMaterialByElement(); }, py::return_value_policy::reference, py::keep_alive<0, 1>()) .def("reassignMaterial", &SolidMechanicsModel::reassignMaterial) .def( "registerNewMaterial", [](SolidMechanicsModel & self, const ID & mat_name, const ID & mat_type, const ID & opt_param) -> decltype(auto) { return self.registerNewMaterial(mat_name, mat_type, opt_param); }, py::arg("material_name"), py::arg("material_type"), py::arg("option") = "", py::return_value_policy::reference) .def("initMaterials", &SolidMechanicsModel::initMaterials) .def("flattenInternal", &SolidMechanicsModel::flattenInternal, py::return_value_policy::reference) .def("inflateInternal", &SolidMechanicsModel::inflateInternal, py::return_value_policy::reference) .def("synchronizeField", &SolidMechanicsModel::synchronizeField, py::arg("synchronization_tag")); } } // namespace akantu diff --git a/python/py_solver.cc b/python/py_solver.cc index 505aaa49b..3a146ede8 100644 --- a/python/py_solver.cc +++ b/python/py_solver.cc @@ -1,144 +1,145 @@ /** * @file py_solver.cc * * @author Nicolas Richart * * @date creation: Tue Sep 29 2020 * @date last modification: Sat Mar 06 2021 * * @brief pybind11 interface to Solver and SparseMatrix * * * @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 "py_solver.hh" #include "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ void register_solvers(py::module & mod) { py::class_(mod, "SparseMatrix") .def("getMatrixType", &SparseMatrix::getMatrixType) .def("size", &SparseMatrix::size) .def("zero", &SparseMatrix::zero) .def("mul", &SparseMatrix::mul) .def("saveProfile", &SparseMatrix::saveProfile) .def("saveMatrix", &SparseMatrix::saveMatrix) .def( "add", [](SparseMatrix & self, UInt i, UInt j) { self.add(i, j); }, "Add entry in the profile") .def( "add", [](SparseMatrix & self, UInt i, UInt j, Real value) { self.add(i, j, value); }, "Add the value to the matrix") .def( "add", [](SparseMatrix & self, SparseMatrix & A, Real alpha) { self.add(A, alpha); }, "Add a matrix to the matrix", py::arg("A"), py::arg("alpha") = 1.) .def("isFinite", &SparseMatrix::isFinite) .def("getRelease", [](const SparseMatrix & self) -> UInt { return self.getRelease(); }) .def("__call__", [](const SparseMatrix & self, UInt i, UInt j) { return self(i, j); }) .def("getRelease", &SparseMatrix::getRelease) .def("applyBoundary", &SparseMatrix::applyBoundary, py::arg("block_val") = 1.); py::class_(mod, "SparseMatrixAIJ") .def("getIRN", &SparseMatrixAIJ::getIRN) .def("getJCN", &SparseMatrixAIJ::getJCN) .def("getA", &SparseMatrixAIJ::getA) .def("copyContent", &SparseMatrixAIJ::copyContent) .def("copyProfile", &SparseMatrixAIJ::copyProfile) .def( "matVecMul", [](const SparseMatrixAIJ & self, const Array & x, Array & y, Real alpha, Real beta) { self.matVecMul(x, y, alpha, beta); }, "Multiply a matrix by a normal vector", py::arg("x"), py::arg("y"), py::arg("alpha") = 1., py::arg("beta") = 0.); py::class_(mod, "SolverVector") .def("zero", &SolverVector::zero) .def( "getValues", [](SolverVector & self) -> decltype(auto) { return static_cast &>(self); }, py::return_value_policy::reference_internal, "Transform this into a vector, Is not copied.") .def("isDistributed", [](const SolverVector & self) { return self.isDistributed(); }) // had to use DUNder methods, otherwise compilation fails .def("__iadd__", &SolverVector::operator+=, py::return_value_policy::reference) .def("__isub__", &SolverVector::operator-=, py::return_value_policy::reference) .def("__imul__", &SolverVector::operator*=, py::return_value_policy::reference) .def("dot", &SolverVector::dot) .def("set", &SolverVector::set) .def("copy", &SolverVector::copy) .def("add", &SolverVector::add) - .def("norm", &SolverVector::norm); + .def("norm", &SolverVector::norm) + .def("normFreeDOFs", &SolverVector::normFreeDOFs); py::class_(mod, "TermToAssemble") .def(py::init()) .def(py::self += Real()) .def_property_readonly("i", &TermsToAssemble::TermToAssemble::i) .def_property_readonly("j", &TermsToAssemble::TermToAssemble::j); py::class_(mod, "TermsToAssemble") .def(py::init()) .def("getDOFIdM", &TermsToAssemble::getDOFIdM) .def("getDOFIdN", &TermsToAssemble::getDOFIdN) .def( "__call__", [](TermsToAssemble & self, UInt i, UInt j, Real val) { auto & term = self(i, j); term = val; return term; }, py::arg("i"), py::arg("j"), py::arg("val") = 0., py::return_value_policy::reference); } } // namespace akantu diff --git a/src/fe_engine/shape_cohesive_inline_impl.hh b/src/fe_engine/shape_cohesive_inline_impl.hh index cf735c74b..3181cadd5 100644 --- a/src/fe_engine/shape_cohesive_inline_impl.hh +++ b/src/fe_engine/shape_cohesive_inline_impl.hh @@ -1,610 +1,610 @@ /** * @file shape_cohesive_inline_impl.hh * * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Fri Feb 03 2012 * @date last modification: Tue Sep 29 2020 * * @brief ShapeCohesive inline implementation * * * @section LICENSE * * Copyright (©) 2010-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 "mesh_iterators.hh" #include "shape_cohesive.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SHAPE_COHESIVE_INLINE_IMPL_HH_ #define AKANTU_SHAPE_COHESIVE_INLINE_IMPL_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ inline ShapeLagrange<_ek_cohesive>::ShapeLagrange(const Mesh & mesh, UInt spatial_dimension, const ID & id) : ShapeLagrangeBase(mesh, spatial_dimension, _ek_cohesive, id) {} #define INIT_SHAPE_FUNCTIONS(type) \ setIntegrationPointsByType(integration_points, ghost_type); \ precomputeShapesOnIntegrationPoints(nodes, ghost_type); \ precomputeShapeDerivativesOnIntegrationPoints(nodes, ghost_type); /* -------------------------------------------------------------------------- */ inline void ShapeLagrange<_ek_cohesive>::initShapeFunctions( const Array & nodes, const Matrix & integration_points, ElementType type, GhostType ghost_type) { AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(INIT_SHAPE_FUNCTIONS); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_cohesive>::computeShapeDerivativesOnIntegrationPoints( const Array & nodes, const Matrix & integration_points, Array & shape_derivatives, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); this->computeShapeDerivativesOnIntegrationPointsLowerDimension( nodes, integration_points, shape_derivatives, ghost_type, filter_elements); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline void ShapeLagrange<_ek_cohesive>::computeShapeDerivativesOnIntegrationPoints( const Array & nodes, const Matrix & integration_points, Array & shape_derivatives, ElementType type, GhostType ghost_type, const Array & filter_elements) const { #define AKANTU_COMPUTE_SHAPES(type) \ computeShapeDerivativesOnIntegrationPoints( \ nodes, integration_points, shape_derivatives, ghost_type, \ filter_elements); AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(AKANTU_COMPUTE_SHAPES); #undef AKANTU_COMPUTE_SHAPES } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_cohesive>::precomputeShapesOnIntegrationPoints( const Array & nodes, GhostType ghost_type) { AKANTU_DEBUG_IN(); InterpolationType itp_type = ElementClassProperty::interpolation_type; Matrix & natural_coords = integration_points(type, ghost_type); UInt size_of_shapes = ElementClass::getShapeSize(); Array & shapes_tmp = shapes.alloc(0, size_of_shapes, itp_type, ghost_type); this->computeShapesOnIntegrationPoints(nodes, natural_coords, shapes_tmp, ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_cohesive>::precomputeShapeDerivativesOnIntegrationPoints( const Array & nodes, GhostType ghost_type) { AKANTU_DEBUG_IN(); InterpolationType itp_type = ElementClassProperty::interpolation_type; Matrix & natural_coords = integration_points(type, ghost_type); UInt size_of_shapesd = ElementClass::getShapeDerivativesSize(); UInt spatial_dimension = mesh.getSpatialDimension(); Array & shapes_derivatives_tmp = shapes_derivatives.alloc(0, size_of_shapesd, itp_type, ghost_type); this->computeShapeDerivativesOnIntegrationPointsLowerDimension( nodes, natural_coords, shapes_derivatives_tmp, ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_cohesive>:: computeShapeDerivativesOnIntegrationPointsLowerDimension( const Array & nodes, const Matrix & integration_points, Array & shape_derivatives, const GhostType & ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt natural_dimension = ElementClass::getNaturalSpaceDimension(); const auto facet_type = Mesh::getFacetType(type); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); UInt nb_points = integration_points.cols(); UInt nb_element = mesh.getConnectivity(type, ghost_type).size(); UInt size_of_shapesd = ElementClass::getShapeDerivativesSize(); AKANTU_DEBUG_ASSERT(shape_derivatives.getNbComponent() == size_of_shapesd, "The shapes_derivatives array does not have the correct " << "number of component"); shape_derivatives.resize(nb_element * nb_points); Array x_el(0, spatial_dimension * nb_nodes_per_element); this->extractNodalToElementField( nodes, x_el, ghost_type, filter_elements); Real * shapesd_val = shape_derivatives.storage(); Array::matrix_iterator x_it = x_el.begin(spatial_dimension, nb_nodes_per_element); Array x_el_equiv(nb_element, natural_dimension * nb_nodes_per_element, 0.); if (spatial_dimension == 3) { AKANTU_DEBUG_ASSERT( facet_type == _triangle_3, "The shape derivatives are calculated only for first order cohesives"); for (auto && data : zip(make_view(x_el, spatial_dimension, nb_nodes_per_element), make_view(x_el_equiv, natural_dimension, nb_nodes_per_element))) { const Matrix & x = std::get<0>(data); auto & x_eq = std::get<1>(data); // compute triangle sides Vector AB = Vector(x(1)) - Vector(x(0)); Vector AC = Vector(x(2)) - Vector(x(0)); Vector BC = Vector(x(2)) - Vector(x(1)); Real a = AB.norm(); Real b = BC.norm(); Real c = AC.norm(); // assign projected coordinates to the nodes of equivalent element // A: [0,0]; B: [a, 0]; C: [xc, yc]; // https://math.stackexchange.com/questions/50227/how-do-i-map-a-3d-triangle-into-2d x_eq(0, 1) = a; auto yc = sqrt((a + b - c) * (a - b + c) * (-a + b + c) * (a + b + c)) / 2 / a; - auto xc = sqrt(c * c - yc * yc); + auto xc = sqrt(std::max(c * c - yc * yc, 0.)); x_eq(0, 2) = xc; x_eq(1, 2) = yc; } } else if (spatial_dimension == 2) { /// array of equivalent coordinates (works for 1D element in 2D only) AKANTU_DEBUG_ASSERT( facet_type == _segment_2, "The shape derivatives are calculated only for first order cohesives"); for (auto && data : zip(make_view(x_el, spatial_dimension, nb_nodes_per_element), make_view(x_el_equiv, natural_dimension, nb_nodes_per_element))) { const Matrix & x = std::get<0>(data); auto & x_eq = std::get<1>(data); for (auto node_nb : arange(1, nb_nodes_per_element)) { Vector delta_x = Vector(x(node_nb)) - Vector(x(0)); x_eq(natural_dimension - 1, node_nb) = delta_x.norm(); } } } if (filter_elements != empty_filter) nb_element = filter_elements.size(); for (auto && data : enumerate( make_view(x_el_equiv, natural_dimension, nb_nodes_per_element))) { auto & elem = std::get<0>(data); const auto & X = std::get<1>(data); if (filter_elements != empty_filter) shapesd_val = shape_derivatives.storage() + filter_elements(elem) * size_of_shapesd * nb_points; Tensor3 B(shapesd_val, natural_dimension, nb_nodes_per_element, nb_points); computeShapeDerivativesOnCPointsByElement(X, integration_points, B); if (filter_elements == empty_filter) shapesd_val += size_of_shapesd * nb_points; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_cohesive>::computeBtD( const Array & Ds, Array & BtDs, GhostType ghost_type, const Array & filter_elements) const { this->computeExtendedBtD(Ds, BtDs, ghost_type, filter_elements); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_cohesive>::computeExtendedBtD( const Array & Ds, Array & AtBtDs, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); auto itp_type = ElementClassProperty::interpolation_type; const auto & shapes_derivatives = this->shapes_derivatives(itp_type, ghost_type); auto natural_dimension = ElementClass::getNaturalSpaceDimension(); auto nb_nodes_per_element = mesh.getNbNodesPerElement(type); Array shapes_derivatives_filtered(0, shapes_derivatives.getNbComponent()); auto && view = make_view(shapes_derivatives, natural_dimension, nb_nodes_per_element / 2); auto B_it = view.begin(); auto B_end = view.end(); if (filter_elements != empty_filter) { FEEngine::filterElementalData(this->mesh, shapes_derivatives, shapes_derivatives_filtered, type, ghost_type, filter_elements); auto && view = make_view(shapes_derivatives_filtered, natural_dimension, nb_nodes_per_element / 2); B_it = view.begin(); B_end = view.end(); } auto A = ExtendingOperators::getAveragingOperator(type); Matrix B_A(natural_dimension, nb_nodes_per_element); for (auto && values : zip(range(B_it, B_end), make_view(Ds, Ds.getNbComponent() / natural_dimension, natural_dimension), make_view(AtBtDs, AtBtDs.getNbComponent() / nb_nodes_per_element, nb_nodes_per_element))) { const auto & B = std::get<0>(values); const auto & D = std::get<1>(values); B_A.mul(B, A); auto & At_Bt_D = std::get<2>(values); // transposed due to the storage layout of B At_Bt_D.template mul(D, B_A); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_cohesive>::computeNtb( const Array & bs, Array & Ntbs, GhostType ghost_type, const Array & filter_elements) const { this->computeExtendedNtb(bs, Ntbs, ghost_type, filter_elements); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_cohesive>::computeExtendedNtb( const Array & bs, Array & CtNtbs, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); CtNtbs.resize(bs.size()); auto itp_type = ElementClassProperty::interpolation_type; auto size_of_shapes = shapes(itp_type, ghost_type).getNbComponent(); auto nb_degree_of_freedom = bs.getNbComponent(); auto nb_nodes_per_element = mesh.getNbNodesPerElement(type); Array shapes_filtered(0, size_of_shapes); auto && view = make_view(shapes(itp_type, ghost_type), 1, size_of_shapes); auto N_it = view.begin(); auto N_end = view.end(); if (filter_elements != empty_filter) { FEEngine::filterElementalData(this->mesh, shapes(itp_type, ghost_type), shapes_filtered, type, ghost_type, filter_elements); auto && view = make_view(shapes_filtered, 1, size_of_shapes); N_it = view.begin(); N_end = view.end(); } auto C = ExtendingOperators::getAveragingOperator(type); Matrix N_C(nb_degree_of_freedom, nb_nodes_per_element); for (auto && values : zip(make_view(bs, nb_degree_of_freedom, 1), range(N_it, N_end), make_view(CtNtbs, nb_degree_of_freedom, nb_nodes_per_element))) { const auto & b = std::get<0>(values); const auto & N = std::get<1>(values); N_C.mul(N, C); auto & CtNtb = std::get<2>(values); CtNtb.template mul(b, N_C); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_cohesive>::extractNodalToElementField( const Array & nodal_f, Array & elemental_f, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); UInt nb_nodes_per_itp_element = ElementClass::getNbNodesPerInterpolationElement(); UInt nb_degree_of_freedom = nodal_f.getNbComponent(); UInt nb_element = this->mesh.getNbElement(type, ghost_type); const auto & conn_array = this->mesh.getConnectivity(type, ghost_type); auto conn = conn_array.begin(conn_array.getNbComponent() / 2, 2); if (filter_elements != empty_filter) { nb_element = filter_elements.size(); } elemental_f.resize(nb_element); Array::matrix_iterator u_it = elemental_f.begin(nb_degree_of_freedom, nb_nodes_per_itp_element); ReduceFunction reduce_function; auto compute = [&](const auto & el) { Matrix & u = *u_it; Matrix el_conn(conn[el]); // compute the average/difference of the nodal field loaded from cohesive // element for (UInt n = 0; n < el_conn.rows(); ++n) { UInt node_plus = el_conn(n, 0); UInt node_minus = el_conn(n, 1); for (UInt d = 0; d < nb_degree_of_freedom; ++d) { Real u_plus = nodal_f(node_plus, d); Real u_minus = nodal_f(node_minus, d); u(d, n) = reduce_function(u_plus, u_minus); } } ++u_it; }; for_each_element(nb_element, filter_elements, compute); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_cohesive>::interpolateOnIntegrationPoints( const Array & in_u, Array & out_uq, UInt nb_degree_of_freedom, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); InterpolationType itp_type = ElementClassProperty::interpolation_type; AKANTU_DEBUG_ASSERT(this->shapes.exists(itp_type, ghost_type), "No shapes for the type " << this->shapes.printType(itp_type, ghost_type)); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); Array u_el(0, nb_degree_of_freedom * nb_nodes_per_element); this->extractNodalToElementField(in_u, u_el, ghost_type, filter_elements); this->template interpolateElementalFieldOnIntegrationPoints( u_el, out_uq, ghost_type, shapes(itp_type, ghost_type), filter_elements); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_cohesive>::variationOnIntegrationPoints( const Array & in_u, Array & nablauq, UInt nb_degree_of_freedom, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); InterpolationType itp_type = ElementClassProperty::interpolation_type; AKANTU_DEBUG_ASSERT( this->shapes_derivatives.exists(itp_type, ghost_type), "No shapes for the type " << this->shapes_derivatives.printType(itp_type, ghost_type)); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); Array u_el(0, nb_degree_of_freedom * nb_nodes_per_element); this->extractNodalToElementField(in_u, u_el, ghost_type, filter_elements); this->template gradientElementalFieldOnIntegrationPoints( u_el, nablauq, ghost_type, shapes_derivatives(itp_type, ghost_type), filter_elements); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_cohesive>::gradientOnIntegrationPoints( const Array & in_u, Array & out_nablauq, UInt nb_degree_of_freedom, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); InterpolationType itp_type = ElementClassProperty::interpolation_type; AKANTU_DEBUG_ASSERT( shapes_derivatives.exists(itp_type, ghost_type), "No shapes derivatives for the type " << shapes_derivatives.printType(itp_type, ghost_type)); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); Array u_el(0, nb_degree_of_freedom * nb_nodes_per_element); this->extractNodalToElementField(in_u, u_el, ghost_type, filter_elements); this->gradientElementalFieldOnIntegrationPoints( u_el, out_nablauq, ghost_type, shapes_derivatives(itp_type, ghost_type), filter_elements); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_cohesive>::computeNormalsOnIntegrationPoints( const Array & u, Array & normals_u, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); UInt nb_element = this->mesh.getNbElement(type, ghost_type); UInt nb_points = this->integration_points(type, ghost_type).cols(); UInt spatial_dimension = this->mesh.getSpatialDimension(); if (filter_elements != empty_filter) { nb_element = filter_elements.size(); } normals_u.resize(nb_points * nb_element); Array tangents_u(0, (spatial_dimension * (spatial_dimension - 1))); if (spatial_dimension > 1) { tangents_u.resize(nb_element * nb_points); this->template variationOnIntegrationPoints( u, tangents_u, spatial_dimension, ghost_type, filter_elements); } Real * tangent = tangents_u.storage(); if (spatial_dimension == 3) { for (auto & normal : make_view(normals_u, spatial_dimension)) { Math::vectorProduct3(tangent, tangent + spatial_dimension, normal.storage()); normal /= normal.norm(); tangent += spatial_dimension * 2; } } else if (spatial_dimension == 2) { for (auto & normal : make_view(normals_u, spatial_dimension)) { Vector a1(tangent, spatial_dimension); normal(0) = -a1(1); normal(1) = a1(0); normal.normalize(); tangent += spatial_dimension; } } else if (spatial_dimension == 1) { const auto facet_type = Mesh::getFacetType(type); const auto & mesh_facets = mesh.getMeshFacets(); const auto & facets = mesh_facets.getSubelementToElement(type, ghost_type); const auto & segments = mesh_facets.getElementToSubelement(facet_type, ghost_type); Real values[2]; for (auto el : arange(nb_element)) { if (filter_elements != empty_filter) { el = filter_elements(el); } for (UInt p = 0; p < 2; ++p) { Element facet = facets(el, p); Element segment = segments(facet.element)[0]; Vector barycenter(values + p, 1); mesh.getBarycenter(segment, barycenter); } Real difference = values[0] - values[1]; AKANTU_DEBUG_ASSERT(difference != 0., "Error in normal computation for cohesive elements"); normals_u(el) = difference / std::abs(difference); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_cohesive>::computeAllignedBasisOnIntegrationPoints( const Array & u, Array & basis_u, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); UInt nb_element = this->mesh.getNbElement(type, ghost_type); UInt nb_points = this->integration_points(type, ghost_type).cols(); UInt spatial_dimension = this->mesh.getSpatialDimension(); if (filter_elements != empty_filter) { nb_element = filter_elements.size(); } basis_u.resize(nb_points * nb_element); Array tangents_u(0, (spatial_dimension * (spatial_dimension - 1))); if (spatial_dimension > 1) { tangents_u.resize(nb_element * nb_points); this->template variationOnIntegrationPoints( u, tangents_u, spatial_dimension, ghost_type, filter_elements); } Real * tangent = tangents_u.storage(); Real * basis_storage = basis_u.storage(); for (auto && data : zip(make_view(basis_u, spatial_dimension, spatial_dimension), make_view(tangents_u, spatial_dimension, (spatial_dimension - 1)))) { auto && basis = std::get<0>(data); auto && tangents = std::get<1>(data); for (auto i : arange(spatial_dimension - 1)) { basis(i) = tangents(i); Vector basis_vec(basis(i), false); basis_vec /= basis_vec.norm(); } if (spatial_dimension == 3) { Vector normal(spatial_dimension); Math::vectorProduct3(tangent, tangent + spatial_dimension, basis_storage + 2 * spatial_dimension); Vector basis_vec(basis(2), false); basis_vec /= basis_vec.norm(); tangent += spatial_dimension * 2; basis_storage += spatial_dimension * 2; } else if (spatial_dimension == 2) { Vector a1(tangent, spatial_dimension); Vector basis_vec(basis(1), false); basis_vec(0) = -a1(1); basis_vec(1) = a1(0); basis_vec /= basis_vec.norm(); tangent += spatial_dimension; } else if (spatial_dimension == 1) { AKANTU_ERROR("1D mesh is currently not supported"); } AKANTU_DEBUG_OUT(); } } } // namespace akantu #endif /* AKANTU_SHAPE_COHESIVE_INLINE_IMPL_HH_ */ diff --git a/src/model/common/non_linear_solver/non_linear_solver_newton_raphson.cc b/src/model/common/non_linear_solver/non_linear_solver_newton_raphson.cc index ca8940469..2463b883d 100644 --- a/src/model/common/non_linear_solver/non_linear_solver_newton_raphson.cc +++ b/src/model/common/non_linear_solver/non_linear_solver_newton_raphson.cc @@ -1,238 +1,243 @@ /** * @file non_linear_solver_newton_raphson.cc * * @author Nicolas Richart * * @date creation: Tue Sep 15 2015 * @date last modification: Tue Mar 30 2021 * * @brief Implementation of the default NonLinearSolver * * * @section LICENSE * * Copyright (©) 2015-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 "non_linear_solver_newton_raphson.hh" #include "communicator.hh" #include "dof_manager_default.hh" #include "solver_callback.hh" #include "solver_vector.hh" #include "sparse_solver_mumps.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ NonLinearSolverNewtonRaphson::NonLinearSolverNewtonRaphson( DOFManagerDefault & dof_manager, const NonLinearSolverType & non_linear_solver_type, const ID & id) : NonLinearSolver(dof_manager, non_linear_solver_type, id), dof_manager(dof_manager), solver(std::make_unique( dof_manager, "J", id + ":sparse_solver")) { this->supported_type.insert(NonLinearSolverType::_newton_raphson_modified); this->supported_type.insert(NonLinearSolverType::_newton_raphson_contact); this->supported_type.insert(NonLinearSolverType::_newton_raphson); this->supported_type.insert(NonLinearSolverType::_linear); this->checkIfTypeIsSupported(); this->registerParam("threshold", convergence_criteria, 1e-10, _pat_parsmod, "Threshold to consider results as converged"); + this->registerParam("threshold_normalized", convergence_criteria_normalized, + _pat_parsmod, + "Threshold to consider results as converged"); this->registerParam("convergence_type", convergence_criteria_type, SolveConvergenceCriteria::_solution, _pat_parsmod, "Type of convergence criteria"); this->registerParam("max_iterations", max_iterations, 10, _pat_parsmod, "Max number of iterations"); this->registerParam("error", error, _pat_readable, "Last reached error"); this->registerParam("nb_iterations", n_iter, _pat_readable, "Last reached number of iterations"); this->registerParam("converged", converged, _pat_readable, "Did last solve converged"); this->registerParam("force_linear_recompute", force_linear_recompute, true, _pat_modifiable, "Force reassembly of the jacobian matrix"); } /* -------------------------------------------------------------------------- */ NonLinearSolverNewtonRaphson::~NonLinearSolverNewtonRaphson() = default; /* ------------------------------------------------------------------------ */ void NonLinearSolverNewtonRaphson::solve(SolverCallback & solver_callback) { solver_callback.beforeSolveStep(); this->dof_manager.updateGlobalBlockedDofs(); solver_callback.predictor(); if (non_linear_solver_type == NonLinearSolverType::_linear and solver_callback.canSplitResidual()) { solver_callback.assembleMatrix("K"); } this->assembleResidual(solver_callback); if (this->non_linear_solver_type == NonLinearSolverType::_newton_raphson_modified || (this->non_linear_solver_type == NonLinearSolverType::_linear && this->force_linear_recompute)) { solver_callback.assembleMatrix("J"); this->force_linear_recompute = false; } this->n_iter = 0; this->converged = false; this->convergence_criteria_normalized = this->convergence_criteria; if (this->convergence_criteria_type == SolveConvergenceCriteria::_residual) { this->converged = this->testConvergence(this->dof_manager.getResidual()); if (this->converged) { return; } auto local_norm = this->normLocalAndMasterDOFs(this->dof_manager.getResidual()); + + local_norm = std::max(local_norm, 1.); this->convergence_criteria_normalized = local_norm * this->convergence_criteria; } do { if (this->non_linear_solver_type == NonLinearSolverType::_newton_raphson or this->non_linear_solver_type == NonLinearSolverType::_newton_raphson_contact) { solver_callback.assembleMatrix("J"); } this->solver->solve(); solver_callback.corrector(); // EventManager::sendEvent(NonLinearSolver::AfterSparseSolve(method)); if (this->convergence_criteria_type == SolveConvergenceCriteria::_residual) { this->assembleResidual(solver_callback); this->converged = this->testConvergence(this->dof_manager.getResidual()); } else { this->converged = this->testConvergence(this->dof_manager.getSolution()); } if (this->convergence_criteria_type == SolveConvergenceCriteria::_solution and not this->converged) { this->assembleResidual(solver_callback); } this->n_iter++; AKANTU_DEBUG_INFO( "[" << this->convergence_criteria_type << "] Convergence iteration " << std::setw(std::log10(this->max_iterations)) << this->n_iter << ": error " << this->error << (this->converged ? " < " : " > ") << this->convergence_criteria); } while (not this->converged and this->n_iter <= this->max_iterations); // this makes sure that you have correct strains and stresses after the // solveStep function (e.g., for dumping) if (this->convergence_criteria_type == SolveConvergenceCriteria::_solution) { this->assembleResidual(solver_callback); } this->converged = this->converged and not(this->n_iter > this->max_iterations); solver_callback.afterSolveStep(this->converged); if (not this->converged) { AKANTU_CUSTOM_EXCEPTION(debug::NLSNotConvergedException( this->convergence_criteria, this->n_iter, this->error)); AKANTU_DEBUG_WARNING("[" << this->convergence_criteria_type << "] Convergence not reached after " << std::setw(std::log10(this->max_iterations)) << this->n_iter << " iteration" << (this->n_iter == 1 ? "" : "s") << "!"); } } /* -------------------------------------------------------------------------- */ bool NonLinearSolverNewtonRaphson::testConvergence( const SolverVector & solver_vector) { AKANTU_DEBUG_IN(); const auto & blocked_dofs = this->dof_manager.getBlockedDOFs(); const Array & array(solver_vector); UInt nb_degree_of_freedoms = array.size(); auto arr_it = array.begin(); auto bld_it = blocked_dofs.begin(); Real norm = 0.; for (UInt n = 0; n < nb_degree_of_freedoms; ++n, ++arr_it, ++bld_it) { bool is_local_node = this->dof_manager.isLocalOrMasterDOF(n); if ((!*bld_it) && is_local_node) { norm += *arr_it * *arr_it; } } dof_manager.getCommunicator().allReduce(norm, SynchronizerOperation::_sum); norm = std::sqrt(norm); AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something went wrong in the solve phase"); this->error = norm; AKANTU_DEBUG_OUT(); return (error < this->convergence_criteria_normalized); } /* -------------------------------------------------------------------------- */ Real NonLinearSolverNewtonRaphson::normLocalAndMasterDOFs( const SolverVector & solver_vector) { AKANTU_DEBUG_IN(); const Array & array(solver_vector); UInt nb_degree_of_freedoms = array.size(); auto arr_it = array.begin(); Real norm = 0.; for (UInt n = 0; n < nb_degree_of_freedoms; ++n, ++arr_it) { bool is_local_node = this->dof_manager.isLocalOrMasterDOF(n); if (is_local_node) { norm += *arr_it * *arr_it; } } dof_manager.getCommunicator().allReduce(norm, SynchronizerOperation::_sum); norm = std::sqrt(norm); AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something went wrong in the solve phase"); AKANTU_DEBUG_OUT(); return norm; } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/heat_transfer/heat_transfer_model.cc b/src/model/heat_transfer/heat_transfer_model.cc index e9e49eb35..c48ccc90c 100644 --- a/src/model/heat_transfer/heat_transfer_model.cc +++ b/src/model/heat_transfer/heat_transfer_model.cc @@ -1,912 +1,912 @@ /** * @file heat_transfer_model.cc * * @author Guillaume Anciaux * @author Lucas Frerot * @author Emil Gallyamov * @author David Simon Kammer * @author Srinivasa Babu Ramisetti * @author Nicolas Richart * @author Rui Wang * * @date creation: Sun May 01 2011 * @date last modification: Fri Apr 09 2021 * * @brief Implementation of HeatTransferModel class * * * @section LICENSE * * Copyright (©) 2010-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 "heat_transfer_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 { namespace heat_transfer { namespace details { class ComputeRhoFunctor { public: ComputeRhoFunctor(const HeatTransferModel & model) : model(model){}; void operator()(Matrix & rho, const Element & el) { auto capacity_it = model.getCapacityArray(el.type, el.ghost_type).begin(); auto density_it = model.getDensityArray(el.type, el.ghost_type).begin(); rho.set(capacity_it[el.element] * density_it[el.element]); } private: const HeatTransferModel & model; }; } // namespace details } // namespace heat_transfer /* -------------------------------------------------------------------------- */ HeatTransferModel::HeatTransferModel(Mesh & mesh, UInt dim, const ID & id, ModelType model_type, std::shared_ptr dof_manager) : Model(mesh, model_type, dof_manager, dim, id), temperature_gradient("temperature_gradient", id), temperature_on_qpoints("temperature_on_qpoints", id), temperature_rate_on_qpoints("temperature_rate_on_qpoints", id), conductivity_on_qpoints("conductivity_on_qpoints", id), k_gradt_on_qpoints("k_gradt_on_qpoints", id), density_array("density_array", id), capacity_array("capacity_array", id), T_ref_array("temperature_reference_array", id), initial_conductivity_array("initial_conductivity_array", id), conductivity_variation_array("conductivity_variation_array", id) { AKANTU_DEBUG_IN(); conductivity = Matrix(this->spatial_dimension, this->spatial_dimension); this->registerDataAccessor(*this); if (this->mesh.isDistributed()) { auto & synchronizer = this->mesh.getElementSynchronizer(); this->registerSynchronizer(synchronizer, SynchronizationTag::_htm_temperature); this->registerSynchronizer(synchronizer, SynchronizationTag::_htm_gradient_temperature); } registerFEEngineObject(id + ":fem", mesh, spatial_dimension); this->mesh.registerDumper("heat_transfer", id, true); this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_regular); this->registerParam("conductivity", conductivity, _pat_parsmod); this->registerParam("conductivity_variation", conductivity_variation, 0., _pat_parsmod); this->registerParam("temperature_reference", T_ref, 0., _pat_parsmod); this->registerParam("capacity", capacity, _pat_parsmod); this->registerParam("density", density, _pat_parsmod); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::initModel() { auto & fem = this->getFEEngine(); fem.initShapeFunctions(_not_ghost); fem.initShapeFunctions(_ghost); temperature_on_qpoints.initialize(fem, _nb_component = 1); temperature_rate_on_qpoints.initialize(fem, _nb_component = 1); temperature_gradient.initialize(fem, _nb_component = spatial_dimension); conductivity_on_qpoints.initialize(fem, _nb_component = spatial_dimension * spatial_dimension); k_gradt_on_qpoints.initialize(fem, _nb_component = spatial_dimension); capacity_array.initialize(fem, _nb_component = 1); density_array.initialize(fem, _nb_component = 1); initial_conductivity_array.initialize(fem, _nb_component = spatial_dimension * spatial_dimension); conductivity_variation_array.initialize(fem, _nb_component = 1); T_ref_array.initialize(fem, _nb_component = 1); this->initBC(*this, *temperature, *increment, *external_heat_rate); } /* -------------------------------------------------------------------------- */ FEEngine & HeatTransferModel::getFEEngineBoundary(const ID & name) { return aka::as_type(getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ HeatTransferModel::~HeatTransferModel() = default; /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleCapacityLumped(GhostType ghost_type) { AKANTU_DEBUG_IN(); auto & fem = getFEEngineClass(); heat_transfer::details::ComputeRhoFunctor compute_rho(*this); for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { fem.assembleFieldLumped(compute_rho, "M", "temperature", this->getDOFManager(), type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MatrixType HeatTransferModel::getMatrixType(const ID & matrix_id) const { if (matrix_id == "K" or matrix_id == "M") { return _symmetric; } return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleConductivityMatrix(); } else if (matrix_id == "M" and need_to_reassemble_capacity) { this->assembleCapacity(); } } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleLumpedMatrix(const ID & matrix_id) { if (matrix_id == "M" and need_to_reassemble_capacity) { this->assembleCapacityLumped(); } } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleResidual() { AKANTU_DEBUG_IN(); this->assembleInternalHeatRate(); this->getDOFManager().assembleToResidual("temperature", *this->external_heat_rate, 1); this->getDOFManager().assembleToResidual("temperature", *this->internal_heat_rate, -1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::predictor() { ++temperature_release; } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleCapacityLumped() { AKANTU_DEBUG_IN(); if (!this->getDOFManager().hasLumpedMatrix("M")) { this->getDOFManager().getNewLumpedMatrix("M"); } this->getDOFManager().zeroLumpedMatrix("M"); assembleCapacityLumped(_not_ghost); assembleCapacityLumped(_ghost); need_to_reassemble_capacity_lumped = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType /*unused*/) { DOFManager & dof_manager = this->getDOFManager(); this->allocNodalField(this->temperature, 1, "temperature"); this->allocNodalField(this->external_heat_rate, 1, "external_heat_rate"); this->allocNodalField(this->internal_heat_rate, 1, "internal_heat_rate"); this->allocNodalField(this->blocked_dofs, 1, "blocked_dofs"); if (!dof_manager.hasDOFs("temperature")) { dof_manager.registerDOFs("temperature", *this->temperature, _dst_nodal); dof_manager.registerBlockedDOFs("temperature", *this->blocked_dofs); } if (time_step_solver_type == TimeStepSolverType::_dynamic || time_step_solver_type == TimeStepSolverType::_dynamic_lumped) { this->allocNodalField(this->temperature_rate, 1, "temperature_rate"); if (!dof_manager.hasDOFsDerivatives("temperature", 1)) { dof_manager.registerDOFsDerivative("temperature", 1, *this->temperature_rate); } } } /* -------------------------------------------------------------------------- */ std::tuple HeatTransferModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_lumped_mass: { return std::make_tuple("explicit_lumped", TimeStepSolverType::_dynamic_lumped); } 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 HeatTransferModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_dynamic_lumped: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["temperature"] = IntegrationSchemeType::_forward_euler; options.solution_type["temperature"] = IntegrationScheme::_temperature_rate; break; } case TimeStepSolverType::_static: { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["temperature"] = IntegrationSchemeType::_pseudo_time; options.solution_type["temperature"] = IntegrationScheme::_not_defined; break; } case TimeStepSolverType::_dynamic: { if (this->method == _explicit_consistent_mass) { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["temperature"] = IntegrationSchemeType::_forward_euler; options.solution_type["temperature"] = IntegrationScheme::_temperature_rate; } else { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["temperature"] = IntegrationSchemeType::_backward_euler; options.solution_type["temperature"] = IntegrationScheme::_temperature; } break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); } return options; } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleConductivityMatrix() { AKANTU_DEBUG_IN(); this->computeConductivityOnQuadPoints(_not_ghost); - if (conductivity_release[_not_ghost] == conductivity_matrix_release) { - return; - } + // if (conductivity_release[_not_ghost] == conductivity_matrix_release) { + // return; + // } AKANTU_DEBUG_ASSERT(this->getDOFManager().hasMatrix("K"), "The K matrix has not been initialized yet."); this->getDOFManager().zeroMatrix("K"); auto & fem = this->getFEEngine(); for (auto && type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { auto nb_element = mesh.getNbElement(type); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_quadrature_points = fem.getNbIntegrationPoints(type); auto bt_d_b = std::make_unique>( nb_element * nb_quadrature_points, nb_nodes_per_element * nb_nodes_per_element, "B^t*D*B"); fem.computeBtDB(conductivity_on_qpoints(type), *bt_d_b, 2, type); /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ auto K_e = std::make_unique>( nb_element, nb_nodes_per_element * nb_nodes_per_element, "K_e"); fem.integrate(*bt_d_b, *K_e, nb_nodes_per_element * nb_nodes_per_element, type); this->getDOFManager().assembleElementalMatricesToMatrix( "K", "temperature", *K_e, type, _not_ghost, _symmetric); } conductivity_matrix_release = conductivity_release[_not_ghost]; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::computeConductivityOnQuadPoints(GhostType ghost_type) { // if already computed once check if need to compute if (not initial_conductivity[ghost_type]) { // if temperature did not change, conductivity will not vary - if (temperature_release == conductivity_release[ghost_type]) { - return; - } - // if conductivity_variation is 0 no need to recompute + // if (temperature_release == conductivity_release[ghost_type]) { + // return; + // } + // // if conductivity_variation is 0 no need to recompute if (conductivity_variation == 0.) { return; } } for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { auto & temperature_interpolated = temperature_on_qpoints(type, ghost_type); // compute the temperature on quadrature points this->getFEEngine().interpolateOnIntegrationPoints( *temperature, temperature_interpolated, 1, type, ghost_type); auto & cond = conductivity_on_qpoints(type, ghost_type); auto & initial_cond = initial_conductivity_array(type, ghost_type); auto & T_ref = T_ref_array(type, ghost_type); auto & cond_var = conductivity_variation_array(type, ghost_type); for (auto && tuple : zip(make_view(cond, spatial_dimension, spatial_dimension), temperature_interpolated, make_view(initial_cond, spatial_dimension, spatial_dimension), T_ref, cond_var)) { auto & C = std::get<0>(tuple); auto & T = std::get<1>(tuple); auto & init_C = std::get<2>(tuple); auto & T_reference = std::get<3>(tuple); auto & C_var = std::get<4>(tuple); C = init_C; Matrix variation(spatial_dimension, spatial_dimension, C_var * (T - T_reference)); C += variation; } } conductivity_release[ghost_type] = temperature_release; initial_conductivity[ghost_type] = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::computeKgradT(GhostType ghost_type) { computeConductivityOnQuadPoints(ghost_type); for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { auto & gradient = temperature_gradient(type, ghost_type); for (auto && values : zip(make_view(conductivity_on_qpoints(type, ghost_type), spatial_dimension, spatial_dimension), make_view(gradient, spatial_dimension), make_view(k_gradt_on_qpoints(type, ghost_type), spatial_dimension))) { const auto & C = std::get<0>(values); const auto & BT = std::get<1>(values); auto & k_BT = std::get<2>(values); k_BT.mul(C, BT); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::computeGradT(GhostType ghost_type) { for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { auto & gradient = temperature_gradient(type, ghost_type); this->getFEEngine().gradientOnIntegrationPoints(*temperature, gradient, 1, type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleInternalHeatRate() { AKANTU_DEBUG_IN(); this->internal_heat_rate->zero(); computeGradT(_not_ghost); // communicate the stresses AKANTU_DEBUG_INFO("Send data for residual assembly"); this->asynchronousSynchronize(SynchronizationTag::_htm_gradient_temperature); assembleInternalHeatRate(_not_ghost); // finalize communications AKANTU_DEBUG_INFO("Wait distant stresses"); this->waitEndSynchronize(SynchronizationTag::_htm_gradient_temperature); assembleInternalHeatRate(_ghost); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleInternalHeatRate(GhostType ghost_type) { AKANTU_DEBUG_IN(); auto & fem = this->getFEEngine(); computeKgradT(ghost_type); for (auto type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto & k_gradt_on_qpoints_vect = k_gradt_on_qpoints(type, ghost_type); UInt nb_quad_points = k_gradt_on_qpoints_vect.size(); Array bt_k_gT(nb_quad_points, nb_nodes_per_element); fem.computeBtD(k_gradt_on_qpoints_vect, bt_k_gT, type, ghost_type); UInt nb_elements = mesh.getNbElement(type, ghost_type); Array int_bt_k_gT(nb_elements, nb_nodes_per_element); fem.integrate(bt_k_gT, int_bt_k_gT, nb_nodes_per_element, type, ghost_type); this->getDOFManager().assembleElementalArrayLocalArray( int_bt_k_gT, *this->internal_heat_rate, type, ghost_type, 1); } } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getStableTimeStep() { AKANTU_DEBUG_IN(); Real el_size; Real min_el_size = std::numeric_limits::max(); Real conductivitymax = conductivity(0, 0); Real densitymin = density; Real capacitymin = capacity; // get maximum conductivity, and minimum density and capacity for (auto type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { for (auto && data : zip(make_view(conductivity_on_qpoints(type), spatial_dimension, spatial_dimension), density_array(type), capacity_array(type))) { auto && conductivity_ip = std::get<0>(data); auto && density_ip = std::get<1>(data); auto && capacity_ip = std::get<2>(data); // get the biggest parameter from k11 until k33// for (UInt i = 0; i < spatial_dimension; i++) { for (UInt j = 0; j < spatial_dimension; j++) { conductivitymax = std::max(conductivity_ip(i, j), conductivitymax); } } // get the smallest density and capacity densitymin = std::min(density_ip, densitymin); capacitymin = std::min(capacity_ip, capacitymin); } } for (auto && type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); Array coord(0, nb_nodes_per_element * spatial_dimension); FEEngine::extractNodalToElementField(mesh, mesh.getNodes(), coord, type, _not_ghost); auto el_coord = coord.begin(spatial_dimension, nb_nodes_per_element); UInt nb_element = mesh.getNbElement(type); for (UInt el = 0; el < nb_element; ++el, ++el_coord) { el_size = getFEEngine().getElementInradius(*el_coord, type); min_el_size = std::min(min_el_size, el_size); } AKANTU_DEBUG_INFO("The minimum element size : " << min_el_size << ", the max conductivity is : " << conductivitymax << ", the min density is : " << densitymin << ", and the min capacity is : " << capacitymin); } Real min_dt = 2. * min_el_size * min_el_size / 4. * densitymin * capacitymin / conductivitymax; mesh.getCommunicator().allReduce(min_dt, SynchronizerOperation::_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ void HeatTransferModel::setTimeStep(Real time_step, const ID & solver_id) { Model::setTimeStep(time_step, solver_id); this->mesh.getDumper("heat_transfer").setTimeStep(time_step); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::readMaterials() { auto sect = this->getParserSection(); if (not std::get<1>(sect)) { const auto & section = std::get<0>(sect); this->parseSection(section); } conductivity_on_qpoints.set(conductivity); initial_conductivity_array.set(conductivity); density_array.set(density); capacity_array.set(capacity); conductivity_variation_array.set(conductivity_variation); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); readMaterials(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assembleCapacity() { AKANTU_DEBUG_IN(); auto ghost_type = _not_ghost; this->getDOFManager().zeroMatrix("M"); auto & fem = getFEEngineClass(); heat_transfer::details::ComputeRhoFunctor rho_functor(*this); for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { fem.assembleFieldMatrix(rho_functor, "M", "temperature", this->getDOFManager(), type, ghost_type); } need_to_reassemble_capacity = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferModel::computeRho(Array & rho, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); FEEngine & fem = this->getFEEngine(); UInt nb_element = mesh.getNbElement(type, ghost_type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); rho.resize(nb_element * nb_quadrature_points); for (auto && data : zip(rho, capacity_array(type, ghost_type), density_array(type, ghost_type))) { auto && rho_ip = std::get<0>(data); auto && capacity_ip = std::get<1>(data); auto && density_ip = std::get<2>(data); rho_ip = capacity_ip * density_ip; } // Real * rho_1_val = rho.storage(); // /// compute @f$ rho @f$ for each nodes of each element // for (UInt el = 0; el < nb_element; ++el) { // for (UInt n = 0; n < nb_quadrature_points; ++n) { // *rho_1_val++ = this->capacity; // } // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::computeThermalEnergyByNode() { AKANTU_DEBUG_IN(); Real ethermal = 0.; for (auto && pair : enumerate(make_view( *internal_heat_rate, internal_heat_rate->getNbComponent()))) { auto n = std::get<0>(pair); auto & heat_rate = std::get<1>(pair); Real heat = 0.; bool is_local_node = mesh.isLocalOrMasterNode(n); bool count_node = is_local_node; for (UInt i = 0; i < heat_rate.size(); ++i) { if (count_node) { heat += heat_rate[i] * this->getTimeStep(); } } ethermal += heat; } mesh.getCommunicator().allReduce(ethermal, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return ethermal; } /* -------------------------------------------------------------------------- */ template void HeatTransferModel::getThermalEnergy( iterator Eth, Array::const_iterator T_it, const Array::const_iterator & T_end, Array::const_iterator capacity_it, Array::const_iterator density_it) const { for (; T_it != T_end; ++T_it, ++Eth, ++capacity_it, ++density_it) { *Eth = *capacity_it * *density_it * *T_it; } } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getThermalEnergy(ElementType type, UInt index) { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type); Vector Eth_on_quarature_points(nb_quadrature_points); auto T_it = this->temperature_on_qpoints(type).begin(); T_it += index * nb_quadrature_points; auto capacity_it = capacity_array(type, _not_ghost).begin(); capacity_it += index * nb_quadrature_points; auto density_it = density_array(type, _not_ghost).begin(); density_it += index * nb_quadrature_points; auto T_end = T_it + nb_quadrature_points; getThermalEnergy(Eth_on_quarature_points.storage(), T_it, T_end, capacity_it, density_it); return getFEEngine().integrate(Eth_on_quarature_points, type, index); } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getThermalEnergy() { Real Eth = 0; auto & fem = getFEEngine(); for (auto && type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { auto nb_element = mesh.getNbElement(type, _not_ghost); auto nb_quadrature_points = fem.getNbIntegrationPoints(type, _not_ghost); Array Eth_per_quad(nb_element * nb_quadrature_points, 1); auto & temperature_interpolated = temperature_on_qpoints(type); // compute the temperature on quadrature points this->getFEEngine().interpolateOnIntegrationPoints( *temperature, temperature_interpolated, 1, type); auto T_it = temperature_interpolated.begin(); auto T_end = temperature_interpolated.end(); auto capacity_it = capacity_array(type, _not_ghost).begin(); auto density_it = density_array(type, _not_ghost).begin(); getThermalEnergy(Eth_per_quad.begin(), T_it, T_end, capacity_it, density_it); Eth += fem.integrate(Eth_per_quad, type); } return Eth; } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getEnergy(const std::string & id) { AKANTU_DEBUG_IN(); Real energy = 0; if (id == "thermal") { energy = getThermalEnergy(); } // reduction sum over all processors mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real HeatTransferModel::getEnergy(const std::string & id, ElementType type, UInt index) { AKANTU_DEBUG_IN(); Real energy = 0.; if (id == "thermal") { energy = getThermalEnergy(type, index); } AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assignPropertyToPhysicalGroup( const std::string & property_name, const std::string & group_name, Real value) { AKANTU_DEBUG_ASSERT(property_name != "conductivity", "Scalar was provided instead of a conductivity matrix"); auto && el_group = mesh.getElementGroup(group_name); auto & fem = this->getFEEngine(); for (auto ghost_type : ghost_types) { for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { auto nb_quadrature_points = fem.getNbIntegrationPoints(type); auto && elements = el_group.getElements(type, ghost_type); Array::scalar_iterator field_it; if (property_name == "density") { field_it = density_array(type, ghost_type).begin(); } else if (property_name == "capacity") { field_it = capacity_array(type, ghost_type).begin(); } else { AKANTU_EXCEPTION(property_name + " is not a valid material property name."); } for (auto && el : elements) { for (auto && qpoint : arange(nb_quadrature_points)) { field_it[el * nb_quadrature_points + qpoint] = value; } } need_to_reassemble_capacity = true; need_to_reassemble_capacity_lumped = true; } } } /* -------------------------------------------------------------------------- */ void HeatTransferModel::assignPropertyToPhysicalGroup( const std::string & property_name, const std::string & group_name, Matrix cond_matrix) { AKANTU_DEBUG_ASSERT(property_name == "conductivity", "When updating material parameters, only conductivity " "accepts matrix as an input"); auto && el_group = mesh.getElementGroup(group_name); auto & fem = this->getFEEngine(); auto dim = this->getMesh().getSpatialDimension(); for (auto ghost_type : ghost_types) { for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { auto nb_quadrature_points = fem.getNbIntegrationPoints(type); auto && elements = el_group.getElements(type, ghost_type); auto init_cond_it = make_view(initial_conductivity_array(type, ghost_type), dim, dim) .begin(); auto cond_on_quad_it = make_view(conductivity_on_qpoints(type, ghost_type), dim, dim) .begin(); for (auto && el : elements) { for (auto && qpoint : arange(nb_quadrature_points)) { init_cond_it[el * nb_quadrature_points + qpoint] = cond_matrix; cond_on_quad_it[el * nb_quadrature_points + qpoint] = cond_matrix; } } } conductivity_release[ghost_type] += 1; } } /* -------------------------------------------------------------------------- */ std::shared_ptr HeatTransferModel::createNodalFieldBool( const std::string & field_name, const std::string & group_name, __attribute__((unused)) bool padding_flag) { std::map *> uint_nodal_fields; uint_nodal_fields["blocked_dofs"] = blocked_dofs.get(); auto field = mesh.createNodalField(uint_nodal_fields[field_name], group_name); return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr HeatTransferModel::createNodalFieldReal( const std::string & field_name, const std::string & group_name, __attribute__((unused)) bool padding_flag) { if (field_name == "capacity_lumped") { AKANTU_EXCEPTION( "Capacity lumped is a nodal field now stored in the DOF manager." "Therefore it cannot be used by a dumper anymore"); } std::map *> real_nodal_fields; real_nodal_fields["temperature"] = temperature.get(); real_nodal_fields["temperature_rate"] = temperature_rate.get(); real_nodal_fields["external_heat_rate"] = external_heat_rate.get(); real_nodal_fields["internal_heat_rate"] = internal_heat_rate.get(); real_nodal_fields["increment"] = increment.get(); std::shared_ptr field = mesh.createNodalField(real_nodal_fields[field_name], group_name); return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr HeatTransferModel::createElementalField( const std::string & field_name, const std::string & group_name, bool /*padding_flag*/, UInt /*spatial_dimension*/, ElementKind element_kind) { std::shared_ptr field; if (field_name == "partitions") { field = mesh.createElementalField( mesh.getConnectivities(), group_name, this->spatial_dimension, element_kind); } else if (field_name == "temperature_gradient") { ElementTypeMap nb_data_per_elem = this->mesh.getNbDataPerElem(temperature_gradient); field = mesh.createElementalField( temperature_gradient, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); } else if (field_name == "conductivity") { ElementTypeMap nb_data_per_elem = this->mesh.getNbDataPerElem(conductivity_on_qpoints); field = mesh.createElementalField( conductivity_on_qpoints, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); } else if (field_name == "density") { ElementTypeMap nb_data_per_elem = this->mesh.getNbDataPerElem(density_array); field = mesh.createElementalField( density_array, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); } else if (field_name == "capacity") { ElementTypeMap nb_data_per_elem = this->mesh.getNbDataPerElem(capacity_array); field = mesh.createElementalField( capacity_array, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); } return field; } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/heat_transfer_interface/heat_transfer_interface_model.cc b/src/model/heat_transfer_interface/heat_transfer_interface_model.cc index a7b42afa2..7edee652f 100644 --- a/src/model/heat_transfer_interface/heat_transfer_interface_model.cc +++ b/src/model/heat_transfer_interface/heat_transfer_interface_model.cc @@ -1,822 +1,900 @@ /** * @file heat_transfer_interface_model.hh * * @author Emil Gallyamov * * @date creation: Thu Apr 13 2023 * @date last modification: Thu Apr 13 2023 * * @brief Implementation of HeatTransferInterfaceModel class * * * @section LICENSE * * Copyright (©) 2010-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 "heat_transfer_interface_model.hh" #include "dumpable_inline_impl.hh" #include "element_class.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_cohesive.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 { namespace heat_transfer_interface { namespace details { class ComputeRhoFunctor { public: ComputeRhoFunctor(const HeatTransferInterfaceModel & model) : model(model){}; void operator()(Matrix & rho, const Element & el) { - auto opening_it = model.getOpening(el.type, el.ghost_type).begin(); + // auto opening_it = model.getOpening(el.type, el.ghost_type).begin(); + // rho.set(model.getCapacityInCrack() * model.getDensityInCrack() * + // opening_it[el.element]); rho.set(model.getCapacityInCrack() * model.getDensityInCrack() * - opening_it[el.element]); + model.getDefaultOpening()); } private: const HeatTransferInterfaceModel & model; }; } // namespace details } // namespace heat_transfer_interface /* -------------------------------------------------------------------------- */ HeatTransferInterfaceModel::HeatTransferInterfaceModel( Mesh & mesh, UInt dim, const ID & id, std::shared_ptr dof_manager) : HeatTransferModel(mesh, dim, id, ModelType::_heat_transfer_interface_model, dof_manager), + temperature_on_qpoints("opening_on_qpoints", id), opening_on_qpoints("opening_on_qpoints", id), - opening_rate("opening_rate", id), k_long_w("k_long_w", id), - k_perp_over_w("k_perp_over_w", id) { + opening_rate("opening_rate", id), k_long_w("k_long_w", id) /*, + k_perp_over_w("k_perp_over_w", id)*/ +{ AKANTU_DEBUG_IN(); registerFEEngineObject("InterfacesFEEngine", mesh, Model::spatial_dimension); this->mesh.registerDumper("heat_interfaces", id); this->mesh.addDumpMeshToDumper("heat_interfaces", mesh, spatial_dimension, _not_ghost, _ek_cohesive); if (this->mesh.isDistributed()) { /// create the distributed synchronizer for cohesive elements this->cohesive_synchronizer = std::make_unique( mesh, "cohesive_distributed_synchronizer"); auto & synchronizer = mesh.getElementSynchronizer(); this->cohesive_synchronizer->split(synchronizer, [](auto && el) { return Mesh::getKind(el.type) == _ek_cohesive; }); this->registerSynchronizer(*cohesive_synchronizer, SynchronizationTag::_htm_gradient_temperature); auto & mesh_facets = this->mesh.getMeshFacets(); auto & facet_synchronizer = mesh_facets.getElementSynchronizer(); const auto & cfacet_synchronizer = facet_synchronizer; // update the cohesive element synchronizer cohesive_synchronizer->updateSchemes([&](auto && scheme, auto && proc, auto && direction) { auto & facet_scheme = cfacet_synchronizer.getCommunications().getScheme(proc, direction); for (auto && facet : facet_scheme) { const auto & cohesive_element = const_cast(mesh_facets) .getElementToSubelement(facet)[1]; if (cohesive_element == ElementNull or cohesive_element.kind() != _ek_cohesive) { continue; } scheme.push_back(cohesive_element); } }); } this->registerParam("longitudinal_conductivity", longitudinal_conductivity, _pat_parsmod); this->registerParam("transversal_conductivity", transversal_conductivity, _pat_parsmod); this->registerParam("default_opening", default_opening, 1.e-5, _pat_parsmod); this->registerParam("capacity_in_crack", capacity_in_crack, _pat_parsmod); this->registerParam("density_in_crack", density_in_crack, _pat_parsmod); this->registerParam("use_opening_rate", use_opening_rate, _pat_parsmod); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferInterfaceModel::initModel() { Parent::initModel(); auto & fem = this->getFEEngine("InterfacesFEEngine"); fem.initShapeFunctions(_not_ghost); fem.initShapeFunctions(_ghost); this->temperature_gradient.initialize(fem, _nb_component = spatial_dimension); this->k_long_w.initialize(fem, _nb_component = (spatial_dimension - 1) * (spatial_dimension - 1)); - this->k_perp_over_w.initialize(fem, _nb_component = 1); + this->temperature_on_qpoints.initialize(fem, _nb_component = 1); this->opening_on_qpoints.initialize(fem, _nb_component = 1); if (use_opening_rate) { opening_rate.initialize(fem, _nb_component = 1); } } /* -------------------------------------------------------------------------- */ HeatTransferInterfaceModel::~HeatTransferInterfaceModel() = default; /* -------------------------------------------------------------------------- */ void HeatTransferInterfaceModel::predictor() { Parent::predictor(); /// to force HeatTransferModel recompute conductivity matrix ++this->conductivity_release[_not_ghost]; /// same for HeatTransferInterfaceModel ++opening_release; } /* -------------------------------------------------------------------------- */ void HeatTransferInterfaceModel::readMaterials() { auto sect = this->getParserSection(); if (not std::get<1>(sect)) { const auto & section = std::get<0>(sect); this->parseSection(section); } this->conductivity_on_qpoints.set(conductivity); this->initial_conductivity_array.set(conductivity); this->density_array.set(density); this->capacity_array.set(capacity); this->conductivity_variation_array.set(conductivity_variation); opening_on_qpoints.set(default_opening); } /* -------------------------------------------------------------------------- */ void HeatTransferInterfaceModel::assembleCapacity() { AKANTU_DEBUG_IN(); auto ghost_type = _not_ghost; Parent::assembleCapacity(); auto & fem_interface = this->getFEEngine("InterfacesFEEngine"); heat_transfer_interface::details::ComputeRhoFunctor rho_functor(*this); for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { fem_interface.assembleFieldMatrix(rho_functor, "M", "temperature", this->getDOFManager(), type, ghost_type); } need_to_reassemble_capacity = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ // void HeatTransferInterfaceModel::assembleCapacityLumped(GhostType ghost_type) // { // AKANTU_DEBUG_IN(); // Parent::assembleCapacityLumped(ghost_type); // auto & fem_interface = // this->getFEEngineClass("InterfacesFEEngine"); // heat_transfer_interface::details::ComputeRhoFunctor compute_rho(*this); // for (auto && type : // mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { // fem_interface.assembleFieldLumped(compute_rho, "M", "temperature", // this->getDOFManager(), type, // ghost_type); // } // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ void HeatTransferInterfaceModel::assembleConductivityMatrix() { AKANTU_DEBUG_IN(); Parent::assembleConductivityMatrix(); this->computeKLongOnQuadPoints(_not_ghost); - this->computeKTransOnQuadPoints(_not_ghost); + // this->computeKTransOnQuadPoints(_not_ghost); - if ((long_conductivity_release[_not_ghost] == - this->crack_conductivity_matrix_release) and - (perp_conductivity_release[_not_ghost] == - this->crack_conductivity_matrix_release)) { - return; - } + // if ((long_conductivity_release[_not_ghost] == + // this->crack_conductivity_matrix_release) and + // (perp_conductivity_release[_not_ghost] == + // this->crack_conductivity_matrix_release)) { + // return; + // } auto & fem_interface = this->getFEEngine("InterfacesFEEngine"); for (auto && type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_cohesive)) { auto nb_element = mesh.getNbElement(type); if (nb_element == 0) continue; auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_quadrature_points = fem_interface.getNbIntegrationPoints(type); auto && shape_derivatives = fem_interface.getShapesDerivatives(type); auto && shapes = fem_interface.getShapes(type); auto A = ExtendingOperators::getAveragingOperator(type); auto C = ExtendingOperators::getDifferencingOperator(type); // tangent_conductivity_matrix.zero(); auto && long_cond = this->k_long_w(type); - auto && perp_cond = this->k_perp_over_w(type); - Array at_bt_k_long_b_a(nb_element * nb_quadrature_points, - nb_nodes_per_element * nb_nodes_per_element, - "A^t*B^t*k_long*B*A"); + // auto && perp_cond = this->k_perp_over_w(type); + Array at_bt_k_long_w_b_a(nb_element * nb_quadrature_points, + nb_nodes_per_element * nb_nodes_per_element, + "A^t*B^t*k_long*B*A"); Array ct_nt_k_perp_n_c(nb_element * nb_quadrature_points, nb_nodes_per_element * nb_nodes_per_element, "C^t*N^t*k_perp*N*C"); Array tangent_conductivity( nb_element * nb_quadrature_points, nb_nodes_per_element * nb_nodes_per_element, "tangent_conductivity"); Matrix B_A(spatial_dimension - 1, nb_nodes_per_element); - Matrix k_long_B_A(spatial_dimension - 1, nb_nodes_per_element); + Matrix k_long_w_B_A(spatial_dimension - 1, nb_nodes_per_element); // Matrix N(spatial_dimension, nb_nodes_per_element); Matrix N_C(1, nb_nodes_per_element); for (auto && data : zip(make_view(long_cond, spatial_dimension - 1, spatial_dimension - 1), make_view(shape_derivatives, spatial_dimension - 1, nb_nodes_per_element / 2), - make_view(at_bt_k_long_b_a, nb_nodes_per_element, + make_view(at_bt_k_long_w_b_a, nb_nodes_per_element, nb_nodes_per_element), - perp_cond, make_view(shapes, 1, nb_nodes_per_element / 2), + /*perp_cond,*/ make_view(shapes, 1, nb_nodes_per_element / 2), make_view(ct_nt_k_perp_n_c, nb_nodes_per_element, nb_nodes_per_element), make_view(tangent_conductivity, nb_nodes_per_element, nb_nodes_per_element))) { // assemble conductivity contribution along crack - const auto & k_long = std::get<0>(data); + const auto & _k_long_w = std::get<0>(data); const auto & B = std::get<1>(data); - auto & At_Bt_k_long_B_A = std::get<2>(data); + auto & At_Bt_k_long_w_B_A = std::get<2>(data); B_A.mul(B, A); - k_long_B_A.mul(k_long, B_A); - At_Bt_k_long_B_A.mul(k_long_B_A, B_A); + k_long_w_B_A.mul(_k_long_w, B_A); + At_Bt_k_long_w_B_A.mul(k_long_w_B_A, B_A); // assemble conductivity contribution perpendicular to the crack - const auto & k_perp = std::get<3>(data); - const auto & N = std::get<4>(data); - auto & Ct_Nt_k_perp_N_C = std::get<5>(data); - auto & tangent = std::get<6>(data); + // const auto & k_perp = std::get<3>(data); + const auto & N = std::get<3>(data); + auto & Ct_Nt_k_perp_N_C = std::get<4>(data); + auto & tangent = std::get<5>(data); N_C.mul(N, C); - Ct_Nt_k_perp_N_C.mul(N_C, N_C, k_perp); + // Ct_Nt_k_perp_N_C.mul(N_C, N_C, k_perp); + Ct_Nt_k_perp_N_C.mul(N_C, N_C, + this->transversal_conductivity); // summing two contributions of tangent operator - tangent = At_Bt_k_long_B_A + Ct_Nt_k_perp_N_C; + tangent = At_Bt_k_long_w_B_A + Ct_Nt_k_perp_N_C; } auto K_e = std::make_unique>( nb_element, nb_nodes_per_element * nb_nodes_per_element, "K_e"); fem_interface.integrate(tangent_conductivity, *K_e, nb_nodes_per_element * nb_nodes_per_element, type); this->getDOFManager().assembleElementalMatricesToMatrix( "K", "temperature", *K_e, type, _not_ghost, _unsymmetric); } this->crack_conductivity_matrix_release = long_conductivity_release[_not_ghost]; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferInterfaceModel::computeGradAndDeltaT(GhostType ghost_type) { auto & fem_interface = this->getFEEngineClass("InterfacesFEEngine"); for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { auto nb_quadrature_points = fem_interface.getNbIntegrationPoints(type); auto nb_element = mesh.getNbElement(type, ghost_type); auto & gradient = temperature_gradient(type, ghost_type); auto surface_gradient = std::make_unique>( nb_element * nb_quadrature_points, spatial_dimension - 1, "surface_gradient"); auto jump = std::make_unique>(nb_element * nb_quadrature_points, 1, "temperature_jump"); this->getFEEngine("InterfacesFEEngine") .gradientOnIntegrationPoints(*temperature, *surface_gradient, 1, type, ghost_type); #define COMPUTE_JUMP(type) \ fem_interface.getShapeFunctions() \ .interpolateOnIntegrationPoints( \ *temperature, *jump, 1, ghost_type); AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_JUMP); #undef COMPUTE_JUMP - for (auto && data : zip(make_view(gradient, spatial_dimension), make_view(*surface_gradient, spatial_dimension - 1), *jump)) { auto & grad = std::get<0>(data); auto & surf_grad = std::get<1>(data); auto & delta = std::get<2>(data); for (auto i : arange(spatial_dimension - 1)) { grad(i) = surf_grad(i); } grad(spatial_dimension - 1) = delta; } } } +/* -------------------------------------------------------------------------- + */ +void HeatTransferInterfaceModel::computeTempOnQpoints(GhostType ghost_type) { + auto & fem_interface = + this->getFEEngineClass("InterfacesFEEngine"); + for (auto && type : + mesh.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { + auto & t_on_qpoints = temperature_on_qpoints(type, ghost_type); + +#define COMPUTE_T(type) \ + fem_interface.getShapeFunctions() \ + .interpolateOnIntegrationPoints( \ + *temperature, t_on_qpoints, 1, ghost_type); + + AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_T); +#undef COMPUTE_JUMP + } +} +/* -------------------------------------------------------------------------- */ +void HeatTransferInterfaceModel::updateNormalOpeningAtQuadraturePoints( + Array positions, GhostType ghost_type) { + auto & fem_interface = + this->getFEEngineClass("InterfacesFEEngine"); + for (auto && type : + mesh.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { + auto nb_quadrature_points = fem_interface.getNbIntegrationPoints(type); + auto nb_element = mesh.getNbElement(type, ghost_type); + + auto & normal_openings = opening_on_qpoints(type, ghost_type); + auto opening_vectors = + std::make_unique>(nb_element * nb_quadrature_points, + spatial_dimension, "opening_vectors"); + auto normal_vectors = std::make_unique>( + nb_element * nb_quadrature_points, spatial_dimension, "normal_vectors"); + +#define COMPUTE_JUMP(type) \ + fem_interface.getShapeFunctions() \ + .interpolateOnIntegrationPoints( \ + positions, *opening_vectors, spatial_dimension, ghost_type); + + AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_JUMP); +#undef COMPUTE_JUMP +#define COMPUTE_NORMALS(type) \ + fem_interface.getShapeFunctions() \ + .computeNormalsOnIntegrationPoints( \ + mesh.getNodes(), *normal_vectors, ghost_type); + AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_NORMALS); +#undef COMPUTE_NORMALS + + for (auto && data : + zip(make_view(*opening_vectors, spatial_dimension), + make_view(*normal_vectors, spatial_dimension), normal_openings)) { + auto & opening = std::get<0>(data); + auto & normal = std::get<1>(data); + auto & normal_opening = std::get<2>(data); + normal_opening = opening.dot(normal); + } + } + // ++opening_release; +} /* -------------------------------------------------------------------------- */ void HeatTransferInterfaceModel::computeKLongOnQuadPoints( GhostType ghost_type) { // if opening did not change, longitudinal conductivity will neither - if (opening_release == long_conductivity_release[ghost_type]) { - return; - } + // if (opening_release == long_conductivity_release[ghost_type]) { + // return; + // } for (auto && type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { auto & opening = opening_on_qpoints(type, ghost_type); auto & long_cond_w = k_long_w(type, ghost_type); UInt dim = spatial_dimension - 1; Matrix long_cond_vect(dim, dim); for (auto i : arange(dim)) { long_cond_vect(i, i) = this->longitudinal_conductivity; } for (auto && tuple : zip(make_view(long_cond_w, dim, dim), opening)) { auto & k = std::get<0>(tuple); - auto & w = std::get<1>(tuple); - + auto w = std::get<1>(tuple); + // limiting the lower limit of w to avoid numeric instability + if (w < this->default_opening) { + w = this->default_opening; + } k = long_cond_vect * w; } } long_conductivity_release[ghost_type] = opening_release; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -void HeatTransferInterfaceModel::computeKTransOnQuadPoints( - GhostType ghost_type) { - - // if opening did not change, longitudinal conductivity will neither - if (opening_release == perp_conductivity_release[ghost_type]) { - return; - } +// void HeatTransferInterfaceModel::computeKTransOnQuadPoints( +// GhostType ghost_type) { - for (auto && type : - mesh.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { - auto & opening = opening_on_qpoints(type, ghost_type); - auto & trans_cond_over_w = k_perp_over_w(type, ghost_type); +// // if opening did not change, longitudinal conductivity will neither +// // if (opening_release == perp_conductivity_release[ghost_type]) { +// // return; +// // } - for (auto && tuple : zip(trans_cond_over_w, opening)) { - auto & k_perp = std::get<0>(tuple); - auto & w = std::get<1>(tuple); +// for (auto && type : +// mesh.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { +// auto & opening = opening_on_qpoints(type, ghost_type); +// auto & trans_cond_over_w = k_perp_over_w(type, ghost_type); + +// for (auto && tuple : zip(trans_cond_over_w, opening)) { +// auto & k_perp = std::get<0>(tuple); +// auto w = std::get<1>(tuple); +// // limiting the lower limit of w to avoid numeric instability +// if (w < this->default_opening) { +// w = this->default_opening; +// } - k_perp = transversal_conductivity / w; - } - } - perp_conductivity_release[ghost_type] = opening_release; +// k_perp = transversal_conductivity / w; +// } +// } +// perp_conductivity_release[ghost_type] = opening_release; - AKANTU_DEBUG_OUT(); -} +// AKANTU_DEBUG_OUT(); +// } /* -------------------------------------------------------------------------- */ void HeatTransferInterfaceModel::assembleInternalHeatRate() { AKANTU_DEBUG_IN(); Parent::assembleInternalHeatRate(); computeGradAndDeltaT(_not_ghost); // communicate the stresses AKANTU_DEBUG_INFO("Send data for residual assembly"); this->asynchronousSynchronize(SynchronizationTag::_htm_gradient_temperature); computeLongHeatRate(_not_ghost); computeTransHeatRate(_not_ghost); computeInertialHeatRate(_not_ghost); // finalize communications AKANTU_DEBUG_INFO("Wait distant stresses"); this->waitEndSynchronize(SynchronizationTag::_htm_gradient_temperature); computeLongHeatRate(_ghost); computeTransHeatRate(_ghost); computeInertialHeatRate(_ghost); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferInterfaceModel::computeLongHeatRate(GhostType ghost_type) { AKANTU_DEBUG_IN(); computeKLongOnQuadPoints(ghost_type); auto & internal_rate = const_cast &>(this->getInternalHeatRate()); auto & fem = this->getFEEngine("InterfacesFEEngine"); for (auto type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { const auto & shapes_derivatives = fem.getShapesDerivatives(type, ghost_type); auto & gradient = temperature_gradient(type, ghost_type); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); auto nb_element = mesh.getNbElement(type, ghost_type); auto A = ExtendingOperators::getAveragingOperator(type); Matrix B_A(spatial_dimension - 1, nb_nodes_per_element); - auto k_long_w_gradT = - std::make_unique>(nb_element * nb_quadrature_points, - spatial_dimension - 1, "k_long_grad_T"); - - auto bt_k_w_gradT = std::make_unique>( - nb_element * nb_quadrature_points, nb_nodes_per_element); - + // auto k_long_w_gradT = + // std::make_unique>(nb_element * nb_quadrature_points, + // spatial_dimension - 1, + // "k_long_grad_T"); + Array k_long_w_gradT(nb_element * nb_quadrature_points, + spatial_dimension - 1); + // auto bt_k_w_gradT = std::make_unique>( + // nb_element * nb_quadrature_points, nb_nodes_per_element); + Array bt_k_w_gradT(nb_element * nb_quadrature_points, + nb_nodes_per_element); for (auto && values : zip(make_view(k_long_w(type, ghost_type), spatial_dimension - 1, spatial_dimension - 1), make_view(gradient, spatial_dimension), - make_view(*k_long_w_gradT, spatial_dimension - 1), + make_view(k_long_w_gradT, spatial_dimension - 1), make_view(shapes_derivatives, spatial_dimension - 1, nb_nodes_per_element / 2), - make_view(*bt_k_w_gradT, nb_nodes_per_element))) { + make_view(bt_k_w_gradT, nb_nodes_per_element))) { const auto & k_w = std::get<0>(values); const auto & full_gradT = std::get<1>(values); Vector surf_gradT(spatial_dimension - 1); for (auto i : arange(spatial_dimension - 1)) { surf_gradT(i) = full_gradT(i); } auto & k_w_gradT = std::get<2>(values); k_w_gradT.mul(k_w, surf_gradT); /// compute @f$B_a t_i@f$ const auto & B = std::get<3>(values); auto & At_Bt_vector = std::get<4>(values); B_A.mul(B, A); At_Bt_vector.template mul(B_A, k_w_gradT); } - auto long_heat_rate = std::make_unique>( - nb_element, nb_nodes_per_element, "long_heat_rate"); + // auto long_heat_rate = std::make_unique>( + // nb_element, nb_nodes_per_element, "long_heat_rate"); + Array long_heat_rate(nb_element, nb_nodes_per_element); - fem.integrate(*bt_k_w_gradT, *long_heat_rate, nb_nodes_per_element, type, + fem.integrate(bt_k_w_gradT, long_heat_rate, nb_nodes_per_element, type, ghost_type); /// assemble this->getDOFManager().assembleElementalArrayLocalArray( - *long_heat_rate, internal_rate, type, ghost_type, 1); + long_heat_rate, internal_rate, type, ghost_type, 1); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void HeatTransferInterfaceModel::computeTransHeatRate(GhostType ghost_type) { AKANTU_DEBUG_IN(); - computeKTransOnQuadPoints(ghost_type); + // computeKTransOnQuadPoints(ghost_type); auto & internal_rate = const_cast &>(this->getInternalHeatRate()); auto & fem = this->getFEEngine("InterfacesFEEngine"); for (auto type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { const auto & shapes = fem.getShapes(type, ghost_type); auto size_of_shapes = shapes.getNbComponent(); auto & gradient = temperature_gradient(type, ghost_type); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); auto nb_element = mesh.getNbElement(type, ghost_type); auto natural_dimension = spatial_dimension - 1; // difference computing operator auto C = ExtendingOperators::getDifferencingOperator(type); auto nt_k_deltaT = std::make_unique>( nb_element * nb_quadrature_points, nb_nodes_per_element); - auto k_perp_deltaT_over_w = std::make_unique>( + auto k_perp_deltaT = std::make_unique>( nb_element * nb_quadrature_points, 1, "k_perp_deltaT"); for (auto && values : - zip(*k_perp_deltaT_over_w, k_perp_over_w(type, ghost_type), + zip(*k_perp_deltaT, /*k_perp_over_w(type, ghost_type),*/ make_view(gradient, spatial_dimension), make_view(shapes, size_of_shapes), make_view(*nt_k_deltaT, nb_nodes_per_element))) { - auto & _k_perp_deltaT_over_w = std::get<0>(values); - auto & _k_perp_over_w = std::get<1>(values); - const auto & full_gradT = std::get<2>(values); - const auto & N = std::get<3>(values); - auto & Nt_k_deltaT = std::get<4>(values); + auto & _k_perp_deltaT = std::get<0>(values); + // auto & _k_perp_over_w = std::get<1>(values); + const auto & full_gradT = std::get<1>(values); + const auto & N = std::get<2>(values); + auto & Nt_k_deltaT = std::get<3>(values); - _k_perp_deltaT_over_w = - _k_perp_over_w * full_gradT(spatial_dimension - 1); + // _k_perp_deltaT_over_w = + // _k_perp_over_w * full_gradT(spatial_dimension - 1); + _k_perp_deltaT = + this->transversal_conductivity * full_gradT(spatial_dimension - 1); - Nt_k_deltaT.mul(C, N, _k_perp_deltaT_over_w); + Nt_k_deltaT.mul(C, N, _k_perp_deltaT); } auto perp_heat_rate = std::make_unique>( nb_element, nb_nodes_per_element, "perp_heat_rate"); fem.integrate(*nt_k_deltaT, *perp_heat_rate, nb_nodes_per_element, type, ghost_type); /// assemble this->getDOFManager().assembleElementalArrayLocalArray( *perp_heat_rate, internal_rate, type, ghost_type, 1); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------*/ void HeatTransferInterfaceModel::computeInertialHeatRate(GhostType ghost_type) { AKANTU_DEBUG_IN(); if (not this->use_opening_rate) return; auto & internal_rate = const_cast &>(this->getInternalHeatRate()); auto & fem_interface = this->getFEEngineClass("InterfacesFEEngine"); for (auto type : mesh.elementTypes(spatial_dimension, ghost_type, _ek_cohesive)) { const auto & shapes = fem_interface.getShapes(type, ghost_type); UInt size_of_shapes = shapes.getNbComponent(); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_quadrature_points = fem_interface.getNbIntegrationPoints(type, ghost_type); auto nb_element = mesh.getNbElement(type, ghost_type); auto natural_dimension = spatial_dimension - 1; // averaging operator auto A = ExtendingOperators::getAveragingOperator(type); auto && opening_rates = this->opening_rate(type, ghost_type); auto * nt_opening_rate = new Array(nb_element * nb_quadrature_points, nb_nodes_per_element); for (auto && values : zip(make_view(shapes, size_of_shapes), make_view(*nt_opening_rate, nb_nodes_per_element), opening_rates)) { const auto & N = std::get<0>(values); auto & Nt_opening_rate = std::get<1>(values); auto & open_rate = std::get<2>(values); Nt_opening_rate.mul(A, N, open_rate); } auto * inertial_heat_rate = new Array(nb_element, nb_nodes_per_element, "inertial_heat_rate"); fem_interface.integrate(*nt_opening_rate, *inertial_heat_rate, nb_nodes_per_element, type, ghost_type); delete nt_opening_rate; /// assemble this->getDOFManager().assembleElementalArrayLocalArray( *inertial_heat_rate, internal_rate, type, ghost_type, 1); + + delete inertial_heat_rate; } AKANTU_DEBUG_OUT(); } -/* -------------------------------------------------------------------------- - */ -const ShapeLagrange<_ek_cohesive> & -HeatTransferInterfaceModel::getShapeFunctionsCohesive() { - return this->getFEEngineClass("InterfacesFEEngine") - .getShapeFunctions(); -} /* -------------------------------------------------------------------------- */ Real HeatTransferInterfaceModel::getStableTimeStep() { AKANTU_DEBUG_IN(); auto HTM_time_step = Parent::getStableTimeStep(); Real el_size; Real min_el_size = std::numeric_limits::max(); Real max_opening_on_qpoint = std::numeric_limits::min(); for (const auto & type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_cohesive)) { auto nb_nodes_per_element = mesh.getNbNodesPerElement(type); auto & opening = opening_on_qpoints(type, _not_ghost); auto nb_quad_per_element = opening.getNbComponent(); auto mesh_dim = this->mesh.getSpatialDimension(); const auto facet_type = Mesh::getFacetType(type); Array coord(0, nb_nodes_per_element * mesh_dim); FEEngine::extractNodalToElementField(mesh, mesh.getNodes(), coord, type, _not_ghost); for (auto && data : zip(make_view(coord, mesh_dim, nb_nodes_per_element), make_view(opening, nb_quad_per_element))) { Matrix & el_coord = std::get<0>(data); Vector & el_opening = std::get<1>(data); el_size = getFEEngine().getElementInradius(el_coord, facet_type); min_el_size = std::min(min_el_size, el_size); max_opening_on_qpoint = std::max(max_opening_on_qpoint, el_opening.norm()); } AKANTU_DEBUG_INFO("The minimum element size : " << min_el_size << " and the max opening is : " << max_opening_on_qpoint); } Real min_dt = 2. * min_el_size * min_el_size / 4 * density_in_crack * capacity_in_crack / longitudinal_conductivity; mesh.getCommunicator().allReduce(min_dt, SynchronizerOperation::_min); min_dt = std::min(min_dt, HTM_time_step); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ void HeatTransferInterfaceModel::setTimeStep(Real time_step, const ID & solver_id) { Parent::setTimeStep(time_step, solver_id); this->mesh.getDumper("heat_interfaces").setTimeStep(time_step); } // /* // -------------------------------------------------------------------------- // */ void HeatTransferInterfaceModel::assignPropertyToPhysicalGroup( // const std::string & property_name, const std::string & group_name, // Real value) { // AKANTU_DEBUG_ASSERT(property_name != "conductivity", // "Scalar was provided instead of a conductivity // matrix"); // auto && el_group = mesh.getElementGroup(group_name); // auto & fem = this->getFEEngine(); // for (auto ghost_type : ghost_types) { // for (auto && type : // mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { // auto nb_quadrature_points = fem.getNbIntegrationPoints(type); // auto && elements = el_group.getElements(type, ghost_type); // Array::scalar_iterator field_it; // if (property_name == "density") { // field_it = density_array(type, ghost_type).begin(); // } else if (property_name == "capacity") { // field_it = capacity_array(type, ghost_type).begin(); // } else { // AKANTU_EXCEPTION(property_name + // " is not a valid material property name."); // } // for (auto && el : elements) { // for (auto && qpoint : arange(nb_quadrature_points)) { // field_it[el * nb_quadrature_points + qpoint] = value; // } // } // need_to_reassemble_capacity = true; // need_to_reassemble_capacity_lumped = true; // } // } // } // /* // -------------------------------------------------------------------------- // */ void HeatTransferInterfaceModel::assignPropertyToPhysicalGroup( // const std::string & property_name, const std::string & group_name, // Matrix cond_matrix) { // AKANTU_DEBUG_ASSERT(property_name == "conductivity", // "When updating material parameters, only conductivity // " "accepts matrix as an input"); // auto && el_group = mesh.getElementGroup(group_name); // auto & fem = this->getFEEngine(); // auto dim = this->getMesh().getSpatialDimension(); // for (auto ghost_type : ghost_types) { // for (auto && type : // mesh.elementTypes(spatial_dimension, ghost_type, _ek_regular)) { // auto nb_quadrature_points = fem.getNbIntegrationPoints(type); // auto && elements = el_group.getElements(type, ghost_type); // auto init_cond_it = // make_view(initial_conductivity_array(type, ghost_type), dim, dim) // .begin(); // auto cond_on_quad_it = // make_view(conductivity_on_qpoints(type, ghost_type), dim, dim) // .begin(); // for (auto && el : elements) { // for (auto && qpoint : arange(nb_quadrature_points)) { // init_cond_it[el * nb_quadrature_points + qpoint] = cond_matrix; // cond_on_quad_it[el * nb_quadrature_points + qpoint] = // cond_matrix; // } // } // } // conductivity_release[ghost_type] += 1; // } // } /* -------------------------------------------------------------------------*/ std::shared_ptr HeatTransferInterfaceModel::createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, UInt spatial_dimension, ElementKind element_kind) { std::shared_ptr field; auto getNbDataPerElem = [&](auto & field) { const auto & fe_engine = getFEEngine("InterfacesFEEngine"); ElementTypeMap res; for (auto ghost_type : ghost_types) { for (auto & type : field.elementTypes(ghost_type)) { auto nb_quadrature_points = fe_engine.getNbIntegrationPoints(type, ghost_type); res(type, ghost_type) = field(type, ghost_type).getNbComponent() * nb_quadrature_points; } } return res; }; if (element_kind == _ek_regular) { field = Parent::createElementalField(field_name, group_name, padding_flag, spatial_dimension, element_kind); } else if (element_kind == _ek_cohesive) { if (field_name == "partitions") { field = mesh.createElementalField( mesh.getConnectivities(), group_name, this->spatial_dimension, element_kind); } else if (field_name == "opening_on_qpoints") { auto nb_data_per_elem = getNbDataPerElem(opening_on_qpoints); field = mesh.createElementalField( opening_on_qpoints, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); + } else if (field_name == "temperature_on_qpoints") { + auto nb_data_per_elem = getNbDataPerElem(temperature_on_qpoints); + field = mesh.createElementalField( + temperature_on_qpoints, group_name, this->spatial_dimension, + element_kind, nb_data_per_elem); } else if (field_name == "temperature_gradient") { auto nb_data_per_elem = getNbDataPerElem(temperature_gradient); field = mesh.createElementalField( temperature_gradient, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); } } return field; } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/heat_transfer_interface/heat_transfer_interface_model.hh b/src/model/heat_transfer_interface/heat_transfer_interface_model.hh index a903d54d0..db2d83bd4 100644 --- a/src/model/heat_transfer_interface/heat_transfer_interface_model.hh +++ b/src/model/heat_transfer_interface/heat_transfer_interface_model.hh @@ -1,245 +1,254 @@ /** * @file heat_transfer_interface_model.hh * * @author Emil Gallyamov * * @date creation: Thu Apr 13 2023 * @date last modification: Thu Apr 13 2023 * * @brief Model of Heat Transfer expanded to heat diffusion along interfaces * * * @section LICENSE * * Copyright (©) 2010-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 "heat_transfer_model.hh" /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_HEAT_TRANSFER_INTERFACE_MODEL_HH_ #define AKANTU_HEAT_TRANSFER_INTERFACE_MODEL_HH_ namespace akantu { class HeatTransferInterfaceModel : public HeatTransferModel { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: using MyFEEngineCohesiveType = FEEngineTemplate; using Parent = HeatTransferModel; HeatTransferInterfaceModel(Mesh & mesh, UInt dim = _all_dimensions, const ID & id = "heat_transfer_interface_model", std::shared_ptr dof_manager = nullptr); ~HeatTransferInterfaceModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// read one material file to instantiate all the materials void readMaterials() override; /// initialize the model void initModel() override; void predictor() override; /* ------------------------------------------------------------------------ */ /* Methods for explicit */ /* ------------------------------------------------------------------------ */ public: /// compute and get the stable time step Real getStableTimeStep() override; /// set the stable timestep void setTimeStep(Real time_step, const ID & solver_id = "") override; public: // /// calculate the lumped capacity vector for heat transfer problem // void assembleCapacityLumped() override; public: /// compute the internal heat flux void assembleInternalHeatRate() override; /// assemble the conductivity matrix void assembleConductivityMatrix() override; // /// compute normals on cohesive elements // void computeNormal(const Array & position, Array & normal, // ElementType type, GhostType ghost_type = _not_ghost); // /// compute basis alliged with cohesive elements // void computeBasis(const Array & position, Array & basis, // ElementType type, GhostType ghost_type = _not_ghost); /// assemble the capacity matrix void assembleCapacity() override; /// compute the capacity on quadrature points void computeRho(Array & rho, ElementType type, GhostType ghost_type); - /// get FEEngine for cohesive elements - const ShapeLagrange<_ek_cohesive> & getShapeFunctionsCohesive(); - public: // /// asign material properties to physical groups // void assignPropertyToPhysicalGroup(const std::string & property_name, // const std::string & group_name, // Real value); // /// asign material properties to physical groups // void assignPropertyToPhysicalGroup(const std::string & property_name, // const std::string & group_name, // Matrix cond_matrix); + /// given nodal positions, this function computes and updates normal openings + /// at quads of cohesives + void updateNormalOpeningAtQuadraturePoints(Array positions, + GhostType ghost_type); + + /// compute temperature on ip for output + void computeTempOnQpoints(GhostType ghost_type); + private: /// compute the integrated longitudinal conductivity matrix (or scalar in /// 2D) for each quadrature point void computeKLongOnQuadPoints(GhostType ghost_type); /// compute the transversal conductivity scalar devided by opening - void computeKTransOnQuadPoints(GhostType ghost_type); - + // void computeKTransOnQuadPoints(GhostType ghost_type); + /// compute temperature gradient along surface and temperature jump on + /// cohesive elements void computeGradAndDeltaT(GhostType ghost_type); /// compute internal heat rate along the crack and assemble it into internal /// forces void computeLongHeatRate(GhostType ghost_type); /// compute internal heat rate perpendicular to the crack and assemble it into /// internal forces void computeTransHeatRate(GhostType ghost_type); /// compute heat rate contributions due to the opening/closure rate /// (dw/dt) void computeInertialHeatRate(GhostType ghost_type); /// compute an averaging operator of size (nb_nodes_per_itp_type x /// nb_nodes_per_cohesive_type) Matrix getAveragingOperator(const ElementType & type, const UInt & nb_degree_of_freedom = 1); /// compute an differencing operator of size (nb_nodes_per_itp_type x /// nb_nodes_per_cohesive_type) Matrix getDifferencingOperator(const ElementType & type, const UInt & nb_degree_of_freedom = 1); // /// compute the thermal energy // Real computeThermalEnergyByNode(); protected: /// calculate the lumped capacity vector for heat transfer problem (w /// ghost type) void assembleCapacityLumped(GhostType ghost_type) override { AKANTU_TO_IMPLEMENT(); }; /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: inline UInt 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; /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: std::shared_ptr createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, UInt spatial_dimension, ElementKind kind) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(TransversalConductivity, transversal_conductivity, Real); AKANTU_GET_MACRO(LongitudinalConductivity, longitudinal_conductivity, Real); AKANTU_GET_MACRO(CapacityInCrack, capacity_in_crack, Real); AKANTU_GET_MACRO(DensityInCrack, density_in_crack, Real); AKANTU_GET_MACRO(DefaultOpening, default_opening, Real); /// get the conductivity on q points - AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(TransversalConductivityOnQpoints, - k_perp_over_w, Real); + // AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(TransversalConductivityOnQpoints, + // k_perp_over_w, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(LongitudinalConductivityOnQpoints, k_long_w, Real); /// get the aperture on q points AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Opening, opening_on_qpoints, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Opening, opening_on_qpoints, Real); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: Real transversal_conductivity; Real longitudinal_conductivity; Real default_opening; Real capacity_in_crack; Real density_in_crack; /// crack opening interpolated on quadrature points ElementTypeMapArray opening_on_qpoints; /// opening rate at quadrature points ElementTypeMapArray opening_rate; + /// opening rate at quadrature points + ElementTypeMapArray temperature_on_qpoints; + /// longitudinal conductivity tensor (or a scalar in 2D) multiplied by opening /// on quadrature points ElementTypeMapArray k_long_w; /// transversal conductivity scalar devided by opening on quadrature points - ElementTypeMapArray k_perp_over_w; + // ElementTypeMapArray k_perp_over_w; /// @brief boolean enabling opening-rate term into internal heat rate bool use_opening_rate{false}; - std::unordered_map long_conductivity_release{{_not_ghost, 0}, - {_ghost, 0}}; + std::unordered_map long_conductivity_release{ + {_not_ghost, UInt(-1)}, {_ghost, UInt(-1)}}; - std::unordered_map perp_conductivity_release{{_not_ghost, 0}, - {_ghost, 0}}; + // std::unordered_map perp_conductivity_release{ + // {_not_ghost, UInt(-1)}, {_ghost, UInt(-1)}}; UInt crack_conductivity_matrix_release{UInt(-1)}; UInt opening_release{0}; /// cohesive elements synchronizer std::unique_ptr cohesive_synchronizer; }; } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "heat_transfer_interface_model_inline_impl.hh" #endif /* AKANTU_HEAT_TRANSFER_INTERFACE_MODEL_HH_ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh index 5918f5ca7..6538154ca 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,604 +1,608 @@ /** * @file solid_mechanics_model.hh * * @author Guillaume Anciaux * @author Daniel Pino Muñoz * @author Nicolas Richart * * @date creation: Tue Jul 27 2010 * @date last modification: Fri Apr 09 2021 * * @brief Model of Solid Mechanics * * * @section LICENSE * * Copyright (©) 2010-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 "boundary_condition.hh" #include "data_accessor.hh" #include "fe_engine.hh" #include "model.hh" #include "non_local_manager_callback.hh" #include "solid_mechanics_model_event_handler.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SOLID_MECHANICS_MODEL_HH_ #define AKANTU_SOLID_MECHANICS_MODEL_HH_ namespace akantu { class Material; class MaterialSelector; class DumperIOHelper; class NonLocalManager; template class IntegratorGauss; template class ShapeLagrange; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ class SolidMechanicsModel : public Model, public DataAccessor, public DataAccessor, public BoundaryCondition, public NonLocalManagerCallback, public EventHandlerManager { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: class NewMaterialElementsEvent : public NewElementsEvent { public: AKANTU_GET_MACRO_NOT_CONST(MaterialList, material, Array &); AKANTU_GET_MACRO(MaterialList, material, const Array &); protected: Array material; }; using MyFEEngineType = FEEngineTemplate; protected: using EventManager = EventHandlerManager; public: SolidMechanicsModel(Mesh & mesh, UInt dim = _all_dimensions, const ID & id = "solid_mechanics_model", std::shared_ptr dof_manager = nullptr, ModelType model_type = ModelType::_solid_mechanics_model); ~SolidMechanicsModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize completely the model void initFullImpl( const ModelOptions & options = SolidMechanicsModelOptions()) override; public: /// initialize all internal arrays for materials virtual void initMaterials(); protected: /// initialize the model void initModel() override; /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; /* ------------------------------------------------------------------------ */ /* Solver interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the stiffness matrix, virtual void assembleStiffnessMatrix(bool need_to_reassemble = false); /// assembles the internal forces in the array internal_forces virtual void assembleInternalForces(); 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 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; /// Callback for the model to instantiate the matricees when needed void initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) override; protected: /* ------------------------------------------------------------------------ */ TimeStepSolverType getDefaultSolverType() const override; /* ------------------------------------------------------------------------ */ ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; public: bool isDefaultSolverExplicit() { return method == _explicit_lumped_mass || method == _explicit_consistent_mass; } protected: /// update the current position vector void updateCurrentPosition(); /* ------------------------------------------------------------------------ */ /* Materials (solid_mechanics_model_material.cc) */ /* ------------------------------------------------------------------------ */ public: /// register an empty material of a given type Material & registerNewMaterial(const ID & mat_name, const ID & mat_type, const ID & opt_param); /// reassigns materials depending on the material selector virtual void reassignMaterial(); /// apply a constant eigen_grad_u on all quadrature points of a given material virtual void applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const ID & material_name, GhostType ghost_type = _not_ghost); protected: /// register a material in the dynamic database Material & registerNewMaterial(const ParserSection & mat_section); /// read the material files to instantiate all the materials void instantiateMaterials(); /// set the element_id_by_material and add the elements to the good materials virtual void assignMaterialToElements(const ElementTypeMapArray * filter = nullptr); /* ------------------------------------------------------------------------ */ /* Mass (solid_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the mass matrix for consistent mass resolutions void assembleMass(); public: /// 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); /// syncronize fields directly by the model (for python wrapper) void synchronizeField(SynchronizationTag tag) { this->synchronize(tag); }; protected: /// fill a vector of rho void computeRho(Array & rho, ElementType type, GhostType ghost_type); /// compute the kinetic energy Real getKineticEnergy(); Real getKineticEnergy(ElementType type, UInt index); /// compute the external work (for impose displacement, the velocity should be /// given too) Real getExternalWork(); /* ------------------------------------------------------------------------ */ /* NonLocalManager inherited members */ /* ------------------------------------------------------------------------ */ protected: void initializeNonLocal() override; void updateDataForNonLocalCriterion(ElementTypeMapReal & criterion) override; void computeNonLocalStresses(GhostType ghost_type) override; void insertIntegrationPointsInNeighborhoods(GhostType ghost_type) override; /// update the values of the non local internal void updateLocalInternal(ElementTypeMapReal & internal_flat, GhostType ghost_type, ElementKind kind) override; /// copy the results of the averaging in the materials void updateNonLocalInternal(ElementTypeMapReal & internal_flat, GhostType ghost_type, ElementKind kind) override; /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: UInt 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; UInt getNbData(const Array & dofs, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) override; protected: void splitElementByMaterial(const Array & elements, std::vector> & elements_per_mat) const; template void splitByMaterial(const Array & elements, Operation && op) const; /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) override; void onNodesRemoved(const Array & element_list, const Array & new_numbering, const RemovedNodesEvent & event) override; void onElementsAdded(const Array & element_list, const NewElementsEvent & event) override; void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) override; void onElementsChanged(const Array & /*unused*/, const Array & /*unused*/, const ElementTypeMapArray & /*unused*/, const ChangedElementsEvent & /*unused*/) override{}; /* ------------------------------------------------------------------------ */ /* Dumpable interface (kept for convenience) and dumper relative functions */ /* ------------------------------------------------------------------------ */ public: virtual void onDump(); //! decide wether a field is a material internal or not bool isInternal(const std::string & field_name, ElementKind element_kind); //! give the amount of data per element virtual ElementTypeMap getInternalDataPerElem(const std::string & field_name, ElementKind kind); //! flatten a given material internal field ElementTypeMapArray & flattenInternal(const std::string & field_name, ElementKind kind, GhostType ghost_type = _not_ghost); //! flatten all the registered material internals void flattenAllRegisteredInternals(ElementKind kind); //! inverse operation of the flatten void inflateInternal(const std::string & field_name, const ElementTypeMapArray & field, ElementKind kind, GhostType ghost_type = _not_ghost); 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, UInt spatial_dimension, ElementKind kind) override; void dump(const std::string & dumper_name) override; void dump(const std::string & dumper_name, UInt step) override; void dump(const std::string & dumper_name, Real time, UInt step) override; void dump() override; void dump(UInt step) override; void dump(Real time, UInt step) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// set the value of the time step void setTimeStep(Real time_step, const ID & solver_id = "") override; /// get the value of the conversion from forces/ mass to acceleration AKANTU_GET_MACRO(F_M2A, f_m2a, Real); /// set the value of the conversion from forces/ mass to acceleration AKANTU_SET_MACRO(F_M2A, f_m2a, Real); /// get the SolidMechanicsModel::displacement array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Displacement, displacement); /// get the SolidMechanicsModel::displacement array AKANTU_GET_MACRO_DEREF_PTR(Displacement, displacement); /// get the SolidMechanicsModel::previous_displacement array AKANTU_GET_MACRO_DEREF_PTR(PreviousDisplacement, previous_displacement); /// get the SolidMechanicsModel::current_position array const Array & getCurrentPosition(); /// get the SolidMechanicsModel::displacement_increment array AKANTU_GET_MACRO_DEREF_PTR(Increment, displacement_increment); /// get the SolidMechanicsModel::displacement_increment array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Increment, displacement_increment); /// get the lumped SolidMechanicsModel::mass array AKANTU_GET_MACRO_DEREF_PTR(Mass, mass); /// get the SolidMechanicsModel::velocity array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Velocity, velocity); /// get the SolidMechanicsModel::velocity array AKANTU_GET_MACRO_DEREF_PTR(Velocity, velocity); /// get the SolidMechanicsModel::acceleration array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(Acceleration, acceleration); /// get the SolidMechanicsModel::acceleration array AKANTU_GET_MACRO_DEREF_PTR(Acceleration, acceleration); /// get the SolidMechanicsModel::external_force array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(ExternalForce, external_force); /// get the SolidMechanicsModel::external_force array AKANTU_GET_MACRO_DEREF_PTR(ExternalForce, external_force); /// get the SolidMechanicsModel::force array (external forces) [[deprecated("Use getExternalForce instead of this function")]] Array & getForce() { return getExternalForce(); } /// get the SolidMechanicsModel::internal_force array (internal forces) AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(InternalForce, internal_force); /// get the SolidMechanicsModel::internal_force array (internal forces) AKANTU_GET_MACRO_DEREF_PTR(InternalForce, internal_force); /// get the SolidMechanicsModel::blocked_dofs array AKANTU_GET_MACRO_DEREF_PTR_NOT_CONST(BlockedDOFs, blocked_dofs); /// get the SolidMechanicsModel::blocked_dofs array AKANTU_GET_MACRO_DEREF_PTR(BlockedDOFs, blocked_dofs); /// get an iterable on the materials inline decltype(auto) getMaterials(); /// get an iterable on the materials inline decltype(auto) getMaterials() const; /// get a particular material (by numerical material index) inline Material & getMaterial(UInt mat_index); /// get a particular material (by numerical material index) inline const Material & getMaterial(UInt mat_index) const; /// get a particular material (by material name) inline Material & getMaterial(const std::string & name); /// get a particular material (by material name) inline const Material & getMaterial(const std::string & name) const; /// get a particular material (by material name) inline const Material & getMaterial(const Element & element) const; /// get a particular material id from is name inline UInt getMaterialIndex(const std::string & name) const; /// give the number of materials inline UInt getNbMaterials() const { return materials.size(); } /// give the material internal index from its id Int getInternalIndexFromID(const ID & id) const; /// compute the stable time step Real getStableTimeStep(); /** * @brief Returns the total energy for a given energy type * * Energy types of SolidMechanicsModel expected as argument are: * - `kinetic` * - `external work` * * Other energy types are passed on to the materials. All materials should * define a `potential` energy type. For additional energy types, see material * documentation. */ Real getEnergy(const std::string & energy_id); /// Compute energy for an element type and material index Real getEnergy(const std::string & energy_id, ElementType type, UInt index); /// Compute energy for an individual element Real getEnergy(const std::string & energy_id, const Element & element) { return getEnergy(energy_id, element.type, element.element); } /// Compute energy for an element group Real getEnergy(const ID & energy_id, const ID & group_id); AKANTU_GET_MACRO(MaterialByElement, material_index, const ElementTypeMapArray &); AKANTU_GET_MACRO(MaterialLocalNumbering, material_local_numbering, const ElementTypeMapArray &); /// vectors containing local material element index for each global element /// index AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialByElement, material_index, UInt); // AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialByElement, material_index, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialLocalNumbering, material_local_numbering, UInt); // AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialLocalNumbering, // material_local_numbering, UInt); AKANTU_GET_MACRO_NOT_CONST(MaterialSelector, material_selector, std::shared_ptr); void setMaterialSelector(std::shared_ptr material_selector) { this->material_selector = std::move(material_selector); } /// Access the non_local_manager interface AKANTU_GET_MACRO(NonLocalManager, *non_local_manager, NonLocalManager &); /// get the FEEngine object to integrate or interpolate on the boundary FEEngine & getFEEngineBoundary(const ID & name = "") override; + AKANTU_GET_MACRO(DisplacementRelease, displacement_release, UInt); + + AKANTU_SET_MACRO(DisplacementRelease, displacement_release, UInt); + protected: /// compute the stable time step Real getStableTimeStep(GhostType ghost_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// release version of the displacement array UInt displacement_release{0}; /// release version of the current_position array UInt current_position_release{0}; /// Check if materials need to recompute the mass array bool need_to_reassemble_lumped_mass{true}; /// Check if materials need to recompute the mass matrix bool need_to_reassemble_mass{true}; /// mapping between material name and material internal id std::map materials_names_to_id; protected: /// conversion coefficient form force/mass to acceleration Real f_m2a{1.0}; /// displacements array std::unique_ptr> displacement; /// displacements array at the previous time step (used in finite deformation) std::unique_ptr> previous_displacement; /// increment of displacement std::unique_ptr> displacement_increment; /// lumped mass array std::unique_ptr> mass; /// velocities array std::unique_ptr> velocity; /// accelerations array std::unique_ptr> acceleration; /// external forces array std::unique_ptr> external_force; /// internal forces array std::unique_ptr> internal_force; /// array specifing if a degree of freedom is blocked or not std::unique_ptr> blocked_dofs; /// array of current position used during update residual std::unique_ptr> current_position; /// Arrays containing the material index for each element ElementTypeMapArray material_index; /// Arrays containing the position in the element filter of the material /// (material's local numbering) ElementTypeMapArray material_local_numbering; /// list of used materials std::vector> materials; /// class defining of to choose a material std::shared_ptr material_selector; using flatten_internal_map = std::map, std::unique_ptr>>; /// tells if the material are instantiated flatten_internal_map registered_internals; /// non local manager std::unique_ptr non_local_manager; /// tells if the material are instantiated bool are_materials_instantiated{false}; friend class Material; template friend class CouplerSolidContactTemplate; }; /* -------------------------------------------------------------------------- */ namespace BC { namespace Neumann { using FromStress = FromHigherDim; using FromTraction = FromSameDim; } // namespace Neumann } // namespace BC } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "parser.hh" #include "solid_mechanics_model_inline_impl.hh" #include "solid_mechanics_model_tmpl.hh" /* -------------------------------------------------------------------------- */ #endif /* AKANTU_SOLID_MECHANICS_MODEL_HH_ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.cc index 7ecd9b191..fae3c374c 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.cc @@ -1,433 +1,433 @@ /** * @file material_cohesive_linear.cc * * @author Mauro Corrado * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Feb 22 2012 * @date last modification: Thu Jan 14 2021 * * @brief Linear irreversible cohesive law of mixed mode loading with * random stress definition for extrinsic type * * * @section LICENSE * * Copyright (©) 2015-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 "material_cohesive_linear.hh" #include "dof_synchronizer.hh" #include "solid_mechanics_model_cohesive.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialCohesiveLinear::MaterialCohesiveLinear( SolidMechanicsModel & model, const ID & id) : MaterialCohesive(model, id), sigma_c_eff("sigma_c_eff", *this), delta_c_eff("delta_c_eff", *this), insertion_stress("insertion_stress", *this) { AKANTU_DEBUG_IN(); this->registerParam("beta", beta, Real(0.), _pat_parsable | _pat_readable, "Beta parameter"); this->registerParam("G_c", G_c, Real(0.), _pat_parsable | _pat_readable, "Mode I fracture energy"); this->registerParam("penalty", penalty, Real(0.), _pat_parsable | _pat_readable, "Penalty coefficient"); this->registerParam("volume_s", volume_s, Real(0.), _pat_parsable | _pat_readable, "Reference volume for sigma_c scaling"); this->registerParam("m_s", m_s, Real(1.), _pat_parsable | _pat_readable, "Weibull exponent for sigma_c scaling"); this->registerParam("kappa", kappa, Real(1.), _pat_parsable | _pat_readable, "Kappa parameter"); this->registerParam( "contact_after_breaking", contact_after_breaking, false, - _pat_parsable | _pat_readable, + _pat_parsable | _pat_modifiable, "Activation of contact when the elements are fully damaged"); this->registerParam("max_quad_stress_insertion", max_quad_stress_insertion, false, _pat_parsable | _pat_readable, "Insertion of cohesive element when stress is high " "enough just on one quadrature point"); this->registerParam("recompute", recompute, false, _pat_parsable | _pat_modifiable, "recompute solution"); this->use_previous_delta_max = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::initMaterial() { AKANTU_DEBUG_IN(); MaterialCohesive::initMaterial(); sigma_c_eff.initialize(1); delta_c_eff.initialize(1); insertion_stress.initialize(spatial_dimension); if (not Math::are_float_equal(delta_c, 0.)) { delta_c_eff.setDefaultValue(delta_c); } else { delta_c_eff.setDefaultValue(2 * G_c / sigma_c); } if (model->getIsExtrinsic()) { scaleInsertionTraction(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::updateInternalParameters() { /// compute scalars beta2_kappa2 = beta * beta / kappa / kappa; beta2_kappa = beta * beta / kappa; if (Math::are_float_equal(beta, 0)) { beta2_inv = 0; } else { beta2_inv = 1. / beta / beta; } } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::scaleInsertionTraction() { AKANTU_DEBUG_IN(); // do nothing if volume_s hasn't been specified by the user if (Math::are_float_equal(volume_s, 0.)) { return; } const Mesh & mesh_facets = model->getMeshFacets(); const auto & fe_engine = model->getFEEngine(); const auto & fe_engine_facet = model->getFEEngine("FacetsFEEngine"); Real base_sigma_c = sigma_c; for (auto && type_facet : mesh_facets.elementTypes(spatial_dimension - 1)) { const Array> & facet_to_element = mesh_facets.getElementToSubelement(type_facet); UInt nb_facet = facet_to_element.size(); UInt nb_quad_per_facet = fe_engine_facet.getNbIntegrationPoints(type_facet); // iterator to modify sigma_c for all the quadrature points of a facet auto sigma_c_iterator = sigma_c(type_facet).begin_reinterpret(nb_quad_per_facet, nb_facet); for (UInt f = 0; f < nb_facet; ++f, ++sigma_c_iterator) { const std::vector & element_list = facet_to_element(f); // compute bounding volume Real volume = 0; auto elem = element_list.begin(); auto elem_end = element_list.end(); for (; elem != elem_end; ++elem) { if (*elem == ElementNull) { continue; } // unit vector for integration in order to obtain the volume UInt nb_quadrature_points = fe_engine.getNbIntegrationPoints(elem->type); Vector unit_vector(nb_quadrature_points, 1); volume += fe_engine.integrate(unit_vector, elem->type, elem->element, elem->ghost_type); } // scale sigma_c *sigma_c_iterator -= base_sigma_c; *sigma_c_iterator *= std::pow(volume_s / volume, 1. / m_s); *sigma_c_iterator += base_sigma_c; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::checkInsertion( bool check_only) { AKANTU_DEBUG_IN(); const Mesh & mesh_facets = model->getMeshFacets(); CohesiveElementInserter & inserter = model->getElementInserter(); for (auto && type_facet : mesh_facets.elementTypes(spatial_dimension - 1)) { ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); const auto & facets_check = inserter.getCheckFacets(type_facet); auto & f_insertion = inserter.getInsertionFacets(type_facet); auto & f_filter = facet_filter(type_facet); auto & sig_c_eff = sigma_c_eff(type_cohesive); auto & del_c = delta_c_eff(type_cohesive); auto & ins_stress = insertion_stress(type_cohesive); auto & trac_old = tractions.previous(type_cohesive); const auto & f_stress = model->getStressOnFacets(type_facet); const auto & sigma_lim = sigma_c(type_facet); UInt nb_quad_facet = model->getFEEngine("FacetsFEEngine").getNbIntegrationPoints(type_facet); #ifndef AKANTU_NDEBUG UInt nb_quad_cohesive = model->getFEEngine("CohesiveFEEngine") .getNbIntegrationPoints(type_cohesive); AKANTU_DEBUG_ASSERT(nb_quad_cohesive == nb_quad_facet, "The cohesive element and the corresponding facet do " "not have the same numbers of integration points"); #endif UInt nb_facet = f_filter.size(); // if (nb_facet == 0) continue; auto sigma_lim_it = sigma_lim.begin(); Matrix stress_tmp(spatial_dimension, spatial_dimension); Matrix normal_traction(spatial_dimension, nb_quad_facet); Vector stress_check(nb_quad_facet); UInt sp2 = spatial_dimension * spatial_dimension; const auto & tangents = model->getTangents(type_facet); const auto & normals = model->getFEEngine("FacetsFEEngine") .getNormalsOnIntegrationPoints(type_facet); auto normal_begin = normals.begin(spatial_dimension); auto tangent_begin = tangents.begin(tangents.getNbComponent()); auto facet_stress_begin = f_stress.begin(spatial_dimension, spatial_dimension * 2); std::vector new_sigmas; std::vector> new_normal_traction; std::vector new_delta_c; // loop over each facet belonging to this material for (UInt f = 0; f < nb_facet; ++f, ++sigma_lim_it) { UInt facet = f_filter(f); // skip facets where check shouldn't be realized if (!facets_check(facet)) { continue; } // compute the effective norm on each quadrature point of the facet for (UInt q = 0; q < nb_quad_facet; ++q) { UInt current_quad = facet * nb_quad_facet + q; const Vector & normal = normal_begin[current_quad]; const Vector & tangent = tangent_begin[current_quad]; const Matrix & facet_stress_it = facet_stress_begin[current_quad]; // compute average stress on the current quadrature point Matrix stress_1(facet_stress_it.storage(), spatial_dimension, spatial_dimension); Matrix stress_2(facet_stress_it.storage() + sp2, spatial_dimension, spatial_dimension); stress_tmp.copy(stress_1); stress_tmp += stress_2; stress_tmp /= 2.; Vector normal_traction_vec(normal_traction(q)); // compute normal and effective stress stress_check(q) = computeEffectiveNorm(stress_tmp, normal, tangent, normal_traction_vec); } // verify if the effective stress overcomes the threshold Real final_stress = stress_check.mean(); if (max_quad_stress_insertion) { final_stress = *std::max_element( stress_check.storage(), stress_check.storage() + nb_quad_facet); } if (final_stress > *sigma_lim_it) { f_insertion(facet) = true; if (check_only) { continue; } // store the new cohesive material parameters for each quadrature // point for (UInt q = 0; q < nb_quad_facet; ++q) { Real new_sigma = stress_check(q); Vector normal_traction_vec(normal_traction(q)); if (spatial_dimension != 3) { normal_traction_vec *= -1.; } new_sigmas.push_back(new_sigma); new_normal_traction.push_back(normal_traction_vec); Real new_delta; // set delta_c in function of G_c or a given delta_c value if (Math::are_float_equal(delta_c, 0.)) { new_delta = 2 * G_c / new_sigma; } else { new_delta = (*sigma_lim_it) / new_sigma * delta_c; } new_delta_c.push_back(new_delta); } } } // update material data for the new elements UInt old_nb_quad_points = sig_c_eff.size(); UInt new_nb_quad_points = new_sigmas.size(); sig_c_eff.resize(old_nb_quad_points + new_nb_quad_points); ins_stress.resize(old_nb_quad_points + new_nb_quad_points); trac_old.resize(old_nb_quad_points + new_nb_quad_points); del_c.resize(old_nb_quad_points + new_nb_quad_points); for (UInt q = 0; q < new_nb_quad_points; ++q) { sig_c_eff(old_nb_quad_points + q) = new_sigmas[q]; del_c(old_nb_quad_points + q) = new_delta_c[q]; for (UInt dim = 0; dim < spatial_dimension; ++dim) { ins_stress(old_nb_quad_points + q, dim) = new_normal_traction[q](dim); trac_old(old_nb_quad_points + q, dim) = new_normal_traction[q](dim); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::computeTraction( const Array & normal, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators auto traction_it = tractions(el_type, ghost_type).begin(spatial_dimension); auto opening_it = opening(el_type, ghost_type).begin(spatial_dimension); auto contact_traction_it = contact_tractions(el_type, ghost_type).begin(spatial_dimension); auto contact_opening_it = contact_opening(el_type, ghost_type).begin(spatial_dimension); auto normal_it = normal.begin(spatial_dimension); auto traction_end = tractions(el_type, ghost_type).end(spatial_dimension); auto sigma_c_it = sigma_c_eff(el_type, ghost_type).begin(); auto delta_max_it = delta_max(el_type, ghost_type).begin(); auto delta_c_it = delta_c_eff(el_type, ghost_type).begin(); auto damage_it = damage(el_type, ghost_type).begin(); auto insertion_stress_it = insertion_stress(el_type, ghost_type).begin(spatial_dimension); Vector normal_opening(spatial_dimension); Vector tangential_opening(spatial_dimension); /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++opening_it, ++normal_it, ++sigma_c_it, ++delta_max_it, ++delta_c_it, ++damage_it, ++contact_traction_it, ++insertion_stress_it, ++contact_opening_it) { Real normal_opening_norm{0}; Real tangential_opening_norm{0}; bool penetration{false}; this->computeTractionOnQuad( *traction_it, *opening_it, *normal_it, *delta_max_it, *delta_c_it, *insertion_stress_it, *sigma_c_it, normal_opening, tangential_opening, normal_opening_norm, tangential_opening_norm, *damage_it, penetration, *contact_traction_it, *contact_opening_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::computeTangentTraction( ElementType el_type, Array & tangent_matrix, const Array & normal, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators auto tangent_it = tangent_matrix.begin(spatial_dimension, spatial_dimension); auto tangent_end = tangent_matrix.end(spatial_dimension, spatial_dimension); auto normal_it = normal.begin(spatial_dimension); auto opening_it = opening(el_type, ghost_type).begin(spatial_dimension); /// NB: delta_max_it points on delta_max_previous, i.e. the /// delta_max related to the solution of the previous incremental /// step auto delta_max_it = delta_max.previous(el_type, ghost_type).begin(); auto sigma_c_it = sigma_c_eff(el_type, ghost_type).begin(); auto delta_c_it = delta_c_eff(el_type, ghost_type).begin(); auto damage_it = damage(el_type, ghost_type).begin(); auto contact_opening_it = contact_opening(el_type, ghost_type).begin(spatial_dimension); Vector normal_opening(spatial_dimension); Vector tangential_opening(spatial_dimension); for (; tangent_it != tangent_end; ++tangent_it, ++normal_it, ++opening_it, ++delta_max_it, ++sigma_c_it, ++delta_c_it, ++damage_it, ++contact_opening_it) { Real normal_opening_norm{0}; Real tangential_opening_norm{0}; bool penetration{false}; this->computeTangentTractionOnQuad( *tangent_it, *delta_max_it, *delta_c_it, *sigma_c_it, *opening_it, *normal_it, normal_opening, tangential_opening, normal_opening_norm, tangential_opening_norm, *damage_it, penetration, *contact_opening_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(cohesive_linear, MaterialCohesiveLinear); } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.cc index 9c5f940be..fafe1e457 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.cc @@ -1,282 +1,280 @@ /** * @file material_cohesive_linear_friction.cc * * @author Mauro Corrado * * @date creation: Tue Jan 12 2016 * @date last modification: Fri Dec 11 2020 * * @brief Linear irreversible cohesive law of mixed mode loading with * random stress definition for extrinsic type * * * @section LICENSE * * Copyright (©) 2015-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 "material_cohesive_linear_friction.hh" #include "solid_mechanics_model_cohesive.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialCohesiveLinearFriction:: MaterialCohesiveLinearFriction(SolidMechanicsModel & model, const ID & id) : MaterialParent(model, id), residual_sliding("residual_sliding", *this), friction_force("friction_force", *this) { AKANTU_DEBUG_IN(); this->registerParam("mu", mu_max, Real(0.), _pat_parsable | _pat_readable, "Maximum value of the friction coefficient"); this->registerParam("penalty_for_friction", friction_penalty, Real(0.), _pat_parsable | _pat_readable, "Penalty parameter for the friction behavior"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearFriction::initMaterial() { AKANTU_DEBUG_IN(); MaterialParent::initMaterial(); friction_force.initialize(spatial_dimension); residual_sliding.initialize(1); residual_sliding.initializeHistory(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearFriction::computeTraction( __attribute__((unused)) const Array & normal, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); residual_sliding.resize(); friction_force.resize(); /// define iterators auto traction_it = this->tractions(el_type, ghost_type).begin(spatial_dimension); auto traction_end = this->tractions(el_type, ghost_type).end(spatial_dimension); auto opening_it = this->opening(el_type, ghost_type).begin(spatial_dimension); - auto previous_opening_it = - this->opening.previous(el_type, ghost_type).begin(spatial_dimension); auto contact_traction_it = this->contact_tractions(el_type, ghost_type).begin(spatial_dimension); auto contact_opening_it = this->contact_opening(el_type, ghost_type).begin(spatial_dimension); auto normal_it = this->normal.begin(spatial_dimension); auto sigma_c_it = this->sigma_c_eff(el_type, ghost_type).begin(); auto delta_max_it = this->delta_max(el_type, ghost_type).begin(); auto delta_max_prev_it = this->delta_max.previous(el_type, ghost_type).begin(); auto delta_c_it = this->delta_c_eff(el_type, ghost_type).begin(); auto damage_it = this->damage(el_type, ghost_type).begin(); auto insertion_stress_it = this->insertion_stress(el_type, ghost_type).begin(spatial_dimension); auto res_sliding_it = this->residual_sliding(el_type, ghost_type).begin(); auto res_sliding_prev_it = this->residual_sliding.previous(el_type, ghost_type).begin(); auto friction_force_it = this->friction_force(el_type, ghost_type).begin(spatial_dimension); Vector normal_opening(spatial_dimension); Vector tangential_opening(spatial_dimension); if (not this->model->isDefaultSolverExplicit()) { this->delta_max(el_type, ghost_type) .copy(this->delta_max.previous(el_type, ghost_type)); } /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++opening_it, ++normal_it, ++sigma_c_it, ++delta_max_it, ++delta_c_it, ++damage_it, ++contact_traction_it, ++insertion_stress_it, ++contact_opening_it, ++delta_max_prev_it, ++res_sliding_it, - ++res_sliding_prev_it, ++friction_force_it, ++previous_opening_it) { + ++res_sliding_prev_it, ++friction_force_it) { Real normal_opening_norm; Real tangential_opening_norm; bool penetration; this->computeTractionOnQuad( *traction_it, *opening_it, *normal_it, *delta_max_it, *delta_c_it, *insertion_stress_it, *sigma_c_it, normal_opening, tangential_opening, normal_opening_norm, tangential_opening_norm, *damage_it, penetration, *contact_traction_it, *contact_opening_it); if (penetration) { /// the friction coefficient mu increases with the damage. It /// equals the maximum value when damage = 1. // Real damage = std::min(*delta_max_prev_it / *delta_c_it, // Real(1.)); Real mu = mu_max; // * damage; /// the definition of tau_max refers to the opening /// (penetration) of the previous incremental step Real tau_max = mu * this->penalty * (std::abs(normal_opening_norm)); Real delta_sliding_norm = std::abs(tangential_opening_norm - *res_sliding_prev_it); /// tau is the norm of the friction force, acting tangentially to the /// surface Real tau = std::min(friction_penalty * delta_sliding_norm, tau_max); if ((tangential_opening_norm - *res_sliding_prev_it) < 0.0) { tau = -tau; } /// from tau get the x and y components of friction, to be added in the /// force vector Vector tangent_unit_vector(spatial_dimension); tangent_unit_vector = tangential_opening / tangential_opening_norm; *friction_force_it = tau * tangent_unit_vector; /// update residual_sliding if (friction_penalty == 0.) { *res_sliding_it = tangential_opening_norm; } else { *res_sliding_it = tangential_opening_norm - (std::abs(tau) / friction_penalty); } } else { friction_force_it->zero(); } *traction_it += *friction_force_it; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearFriction::computeTangentTraction( ElementType el_type, Array & tangent_matrix, __attribute__((unused)) const Array & normal, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators auto tangent_it = tangent_matrix.begin(spatial_dimension, spatial_dimension); auto tangent_end = tangent_matrix.end(spatial_dimension, spatial_dimension); auto normal_it = this->normal.begin(spatial_dimension); auto opening_it = this->opening(el_type, ghost_type).begin(spatial_dimension); auto previous_opening_it = this->opening.previous(el_type, ghost_type).begin(spatial_dimension); /** * NB: delta_max_it points on delta_max_previous, i.e. the * delta_max related to the solution of the previous incremental * step */ auto delta_max_it = this->delta_max.previous(el_type, ghost_type).begin(); auto sigma_c_it = this->sigma_c_eff(el_type, ghost_type).begin(); auto delta_c_it = this->delta_c_eff(el_type, ghost_type).begin(); auto damage_it = this->damage(el_type, ghost_type).begin(); auto contact_opening_it = this->contact_opening(el_type, ghost_type).begin(spatial_dimension); auto res_sliding_prev_it = this->residual_sliding.previous(el_type, ghost_type).begin(); Vector normal_opening(spatial_dimension); Vector tangential_opening(spatial_dimension); for (; tangent_it != tangent_end; ++tangent_it, ++normal_it, ++opening_it, ++previous_opening_it, ++delta_max_it, ++sigma_c_it, ++delta_c_it, ++damage_it, ++contact_opening_it, ++res_sliding_prev_it) { Real normal_opening_norm; Real tangential_opening_norm; bool penetration; this->computeTangentTractionOnQuad( *tangent_it, *delta_max_it, *delta_c_it, *sigma_c_it, *opening_it, *normal_it, normal_opening, tangential_opening, normal_opening_norm, tangential_opening_norm, *damage_it, penetration, *contact_opening_it); if (penetration) { // Real damage = std::min(*delta_max_it / *delta_c_it, Real(1.)); Real mu = mu_max; // * damage; Real normal_opening_prev_norm = std::min(previous_opening_it->dot(*normal_it), Real(0.)); // Vector normal_opening_prev = (*normal_it); // normal_opening_prev *= normal_opening_prev_norm; Real tau_max = mu * this->penalty * (std::abs(normal_opening_prev_norm)); Real delta_sliding_norm = std::abs(tangential_opening_norm - *res_sliding_prev_it); // tau is the norm of the friction force, acting tangentially to the // surface Real tau = std::min(friction_penalty * delta_sliding_norm, tau_max); if (tau < tau_max && tau_max > Math::getTolerance()) { Matrix I(spatial_dimension, spatial_dimension); I.eye(1.); Matrix n_outer_n(spatial_dimension, spatial_dimension); n_outer_n.outerProduct(*normal_it, *normal_it); Matrix nn(n_outer_n); I -= nn; *tangent_it += I * friction_penalty; } } // check if the tangential stiffness matrix is symmetric // for (UInt h = 0; h < spatial_dimension; ++h){ // for (UInt l = h; l < spatial_dimension; ++l){ // if (l > h){ // Real k_ls = (*tangent_it)[spatial_dimension*h+l]; // Real k_us = (*tangent_it)[spatial_dimension*l+h]; // // std::cout << "k_ls = " << k_ls << std::endl; // // std::cout << "k_us = " << k_us << std::endl; // if (std::abs(k_ls) > 1e-13 && std::abs(k_us) > 1e-13){ // Real error = std::abs((k_ls - k_us) / k_us); // if (error > 1e-10){ // std::cout << "non symmetric cohesive matrix" << std::endl; // // std::cout << "error " << error << std::endl; // } // } // } // } // } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(cohesive_linear_friction, MaterialCohesiveLinearFriction); } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_inline_impl.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_inline_impl.hh index 6e8609fe8..f9b6fca8f 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_inline_impl.hh +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_inline_impl.hh @@ -1,269 +1,270 @@ /** * @file material_cohesive_linear_inline_impl.hh * * @author Mauro Corrado * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Wed Apr 22 2015 * @date last modification: Thu Jan 14 2021 * * @brief Inline functions of the MaterialCohesiveLinear * * * @section LICENSE * * Copyright (©) 2015-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 "material_cohesive_linear.hh" #include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MATERIAL_COHESIVE_LINEAR_INLINE_IMPL_HH_ #define AKANTU_MATERIAL_COHESIVE_LINEAR_INLINE_IMPL_HH_ /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template inline Real MaterialCohesiveLinear::computeEffectiveNorm( const Matrix & stress, const Vector & normal, const Vector & tangent, Vector & normal_traction) const { normal_traction.mul(stress, normal); Real normal_contrib = normal_traction.dot(normal); /// in 3D tangential components must be summed Real tangent_contrib = 0; if (dim == 2) { Real tangent_contrib_tmp = normal_traction.dot(tangent); tangent_contrib += tangent_contrib_tmp * tangent_contrib_tmp; } else if (dim == 3) { for (UInt s = 0; s < dim - 1; ++s) { const Vector tangent_v(tangent.storage() + s * dim, dim); Real tangent_contrib_tmp = normal_traction.dot(tangent_v); tangent_contrib += tangent_contrib_tmp * tangent_contrib_tmp; } } tangent_contrib = std::sqrt(tangent_contrib); normal_contrib = std::max(Real(0.), normal_contrib); return std::sqrt(normal_contrib * normal_contrib + tangent_contrib * tangent_contrib * beta2_inv); } /* -------------------------------------------------------------------------- */ template inline void MaterialCohesiveLinear::computeTractionOnQuad( Vector & traction, Vector & opening, const Vector & normal, Real & delta_max, const Real & delta_c, const Vector & insertion_stress, const Real & sigma_c, Vector & normal_opening, Vector & tangential_opening, Real & normal_opening_norm, Real & tangential_opening_norm, Real & damage, bool & penetration, Vector & contact_traction, Vector & contact_opening) { /// compute normal and tangential opening vectors normal_opening_norm = opening.dot(normal); normal_opening = normal; normal_opening *= normal_opening_norm; tangential_opening = opening; tangential_opening -= normal_opening; tangential_opening_norm = tangential_opening.norm(); /** * compute effective opening displacement * @f$ \delta = \sqrt{ * \frac{\beta^2}{\kappa^2} \Delta_t^2 + \Delta_n^2 } @f$ */ Real delta = tangential_opening_norm * tangential_opening_norm * this->beta2_kappa2; penetration = normal_opening_norm / delta_c < -Math::getTolerance(); // penetration = normal_opening_norm < 0.; if (not this->contact_after_breaking and Math::are_float_equal(damage, 1.)) { penetration = false; } if (penetration) { /// use penalty coefficient in case of penetration contact_traction = normal_opening; contact_traction *= this->penalty; contact_opening = normal_opening; /// don't consider penetration contribution for delta opening = tangential_opening; normal_opening.zero(); } else { delta += normal_opening_norm * normal_opening_norm; contact_traction.zero(); contact_opening.zero(); } delta = std::sqrt(delta); /// update maximum displacement and damage delta_max = std::max(delta_max, delta); + delta_max = std::min(delta_max, delta_c); damage = std::min(delta_max / delta_c, Real(1.)); /** * Compute traction @f$ \mathbf{T} = \left( * \frac{\beta^2}{\kappa} \Delta_t \mathbf{t} + \Delta_n * \mathbf{n} \right) \frac{\sigma_c}{\delta} \left( 1- * \frac{\delta}{\delta_c} \right)@f$ */ if (Math::are_float_equal(damage, 1.)) { traction.zero(); } else if (Math::are_float_equal(damage, 0.)) { if (penetration) { traction.zero(); } else { traction = insertion_stress; } } else { traction = tangential_opening; traction *= this->beta2_kappa; traction += normal_opening; AKANTU_DEBUG_ASSERT(delta_max != 0., "Division by zero, tolerance might be too low"); traction *= sigma_c / delta_max * (1. - damage); } } /* -------------------------------------------------------------------------- */ template inline void MaterialCohesiveLinear::computeTangentTractionOnQuad( Matrix & tangent, Real & delta_max, const Real & delta_c, const Real & sigma_c, Vector & opening, const Vector & normal, Vector & normal_opening, Vector & tangential_opening, Real & normal_opening_norm, Real & tangential_opening_norm, Real & damage, bool & penetration, Vector & contact_opening) { /** * During the update of the residual the interpenetrations are * stored in the array "contact_opening", therefore, in the case * of penetration, in the array "opening" there are only the * tangential components. */ opening += contact_opening; /// compute normal and tangential opening vectors normal_opening_norm = opening.dot(normal); normal_opening = normal; normal_opening *= normal_opening_norm; tangential_opening = opening; tangential_opening -= normal_opening; tangential_opening_norm = tangential_opening.norm(); Real delta = tangential_opening_norm * tangential_opening_norm * this->beta2_kappa2; penetration = normal_opening_norm < 0.0; if (not this->contact_after_breaking and Math::are_float_equal(damage, 1.)) { penetration = false; } Real derivative = 0; // derivative = d(t/delta)/ddelta Real t = 0; Matrix n_outer_n(spatial_dimension, spatial_dimension); n_outer_n.outerProduct(normal, normal); if (penetration) { /// stiffness in compression given by the penalty parameter tangent += n_outer_n; tangent *= penalty; opening = tangential_opening; normal_opening_norm = opening.dot(normal); normal_opening = normal; normal_opening *= normal_opening_norm; } else { delta += normal_opening_norm * normal_opening_norm; } delta = std::sqrt(delta); /** * Delta has to be different from 0 to have finite values of * tangential stiffness. At the element insertion, delta = * 0. Therefore, a fictictious value is defined, for the * evaluation of the first value of K. */ if (delta < Math::getTolerance()) { delta = delta_c / 1000.; } if (delta >= delta_max) { if (delta <= delta_c) { derivative = -sigma_c / (delta * delta); t = sigma_c * (1 - delta / delta_c); } else { derivative = 0.; t = 0.; } } else if (delta < delta_max) { Real tmax = sigma_c * (1 - delta_max / delta_c); t = tmax / delta_max * delta; } /// computation of the derivative of the constitutive law (dT/ddelta) Matrix I(spatial_dimension, spatial_dimension); I.eye(this->beta2_kappa); Matrix nn(n_outer_n); nn *= (1. - this->beta2_kappa); nn += I; nn *= t / delta; Vector t_tilde(normal_opening); t_tilde *= (1. - this->beta2_kappa2); Vector mm(opening); mm *= this->beta2_kappa2; t_tilde += mm; Vector t_hat(normal_opening); t_hat += this->beta2_kappa * tangential_opening; Matrix prov(spatial_dimension, spatial_dimension); prov.outerProduct(t_hat, t_tilde); prov *= derivative / delta; prov += nn; Matrix prov_t = prov.transpose(); tangent += prov_t; } /* -------------------------------------------------------------------------- */ } // namespace akantu /* -------------------------------------------------------------------------- */ #endif // AKANTU_MATERIAL_COHESIVE_LINEAR_INLINE_IMPL_HH_ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.hh index 9f927aaca..cb28d6be3 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.hh +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/material_cohesive.hh @@ -1,237 +1,237 @@ /** * @file material_cohesive.hh * * @author Nicolas Richart * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * * @date creation: Fri Jun 18 2010 * @date last modification: Thu Jan 14 2021 * * @brief Specialization of the material class for cohesive elements * * * @section LICENSE * * Copyright (©) 2010-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 "material.hh" /* -------------------------------------------------------------------------- */ #include "cohesive_internal_field.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MATERIAL_COHESIVE_HH_ #define AKANTU_MATERIAL_COHESIVE_HH_ /* -------------------------------------------------------------------------- */ namespace akantu { class SolidMechanicsModelCohesive; } namespace akantu { class MaterialCohesive : public Material { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: using MyFEEngineCohesiveType = FEEngineTemplate; public: MaterialCohesive(SolidMechanicsModel & model, const ID & id = ""); ~MaterialCohesive() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize the material computed parameter void initMaterial() override; /// compute tractions (including normals and openings) void computeTraction(GhostType ghost_type = _not_ghost); /// assemble residual void assembleInternalForces(GhostType ghost_type = _not_ghost) override; /// check stress for cohesive elements' insertion, by default it /// also updates the cohesive elements' data virtual void checkInsertion(bool /*check_only*/ = false) { AKANTU_TO_IMPLEMENT(); } /// interpolate stress on given positions for each element (empty /// implemantation to avoid the generic call to be done on cohesive elements) virtual void interpolateStress(const ElementType /*type*/, Array & /*result*/) {} /// compute the stresses void computeAllStresses(GhostType /*ghost_type*/ = _not_ghost) final{}; // add the facet to be handled by the material UInt addFacet(const Element & element); + /// compute the normal + void computeNormal(const Array & position, Array & normal, + ElementType type, GhostType ghost_type); + protected: virtual void computeTangentTraction(ElementType /*el_type*/, Array & /*tangent_matrix*/, const Array & /*normal*/, GhostType /*ghost_type*/ = _not_ghost) { AKANTU_TO_IMPLEMENT(); } - /// compute the normal - void computeNormal(const Array & position, Array & normal, - ElementType type, GhostType ghost_type); - /// compute the opening void computeOpening(const Array & displacement, Array & opening, ElementType type, GhostType ghost_type); template void computeNormal(const Array & position, Array & normal, GhostType ghost_type); /// assemble stiffness void assembleStiffnessMatrix(GhostType ghost_type) override; /// constitutive law virtual void computeTraction(const Array & normal, ElementType el_type, GhostType ghost_type = _not_ghost) = 0; /// parallelism functions inline UInt 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; protected: void updateEnergies(ElementType el_type) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the opening AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Opening, opening, Real); /// get the traction AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Traction, tractions, Real); /// get damage AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Damage, damage, Real); /// get facet filter AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(FacetFilter, facet_filter, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(FacetFilter, facet_filter, UInt); AKANTU_GET_MACRO(FacetFilter, facet_filter, const ElementTypeMapArray &); AKANTU_GET_MACRO(NormalsAtQuads, normal, const Array &); // AKANTU_GET_MACRO(ElementFilter, element_filter, const // ElementTypeMapArray &); /// compute reversible energy Real getReversibleEnergy(); /// compute dissipated energy Real getDissipatedEnergy(); /// compute contact energy Real getContactEnergy(); /// get energy Real getEnergy(const std::string & type) override; /// return the energy (identified by id) for the provided element Real getEnergy(const std::string & energy_id, ElementType type, UInt index) override { return Material::getEnergy(energy_id, type, index); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// list of facets assigned to this material ElementTypeMapArray facet_filter; /// Link to the cohesive fem object in the model FEEngine & fem_cohesive; private: /// reversible energy by quadrature point CohesiveInternalField reversible_energy; /// total energy by quadrature point CohesiveInternalField total_energy; protected: /// opening in all elements and quadrature points CohesiveInternalField opening; /// traction in all elements and quadrature points CohesiveInternalField tractions; /// traction due to contact CohesiveInternalField contact_tractions; /// normal openings for contact tractions CohesiveInternalField contact_opening; /// maximum displacement CohesiveInternalField delta_max; /// tell if the previous delta_max state is needed (in iterative schemes) bool use_previous_delta_max; /// tell if the previous opening state is needed (in iterative schemes) bool use_previous_opening; /// damage CohesiveInternalField damage; /// pointer to the solid mechanics model for cohesive elements SolidMechanicsModelCohesive * model; /// critical stress RandomInternalField sigma_c; /// critical displacement Real delta_c; /// array to temporarily store the normals Array normal; }; } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "cohesive_internal_field_tmpl.hh" #include "material_cohesive_inline_impl.hh" #endif /* AKANTU_MATERIAL_COHESIVE_HH_ */ diff --git a/src/solver/solver_vector.hh b/src/solver/solver_vector.hh index 237f10ece..0e850d2b6 100644 --- a/src/solver/solver_vector.hh +++ b/src/solver/solver_vector.hh @@ -1,103 +1,104 @@ /** * @file solver_vector.hh * * @author Nicolas Richart * * @date creation: Thu Feb 21 2013 * @date last modification: Tue May 26 2020 * * @brief Solver vector interface base 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 "aka_array.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SOLVER_VECTOR_HH_ #define AKANTU_SOLVER_VECTOR_HH_ namespace akantu { class DOFManager; } namespace akantu { class SolverVector { public: SolverVector(DOFManager & dof_manager, const ID & id = "solver_vector") : id(id), _dof_manager(dof_manager) {} SolverVector(const SolverVector & vector, const ID & id = "solver_vector") : id(id), _dof_manager(vector._dof_manager) {} virtual ~SolverVector() = default; // resize the vector to the size of the problem virtual void resize() = 0; // clear the vector virtual void set(Real val) = 0; void zero() { this->set({}); } virtual operator const Array &() const = 0; virtual Int size() const = 0; virtual Int localSize() const = 0; virtual SolverVector & operator+=(const SolverVector & y) = 0; virtual SolverVector & operator-=(const SolverVector & y) = 0; virtual SolverVector & operator=(const SolverVector & y) = 0; virtual void copy(const SolverVector & y) = 0; virtual SolverVector & operator*=(const Real & alpha) = 0; virtual void add(const SolverVector & y, const Real & alpha) = 0; virtual Real dot(const SolverVector & y) const = 0; /// computes l2 norm of a vector virtual Real norm() const = 0; + virtual Real normFreeDOFs() const = 0; UInt & release() { return release_; } UInt release() const { return release_; } virtual void printself(std::ostream & stream, int indent = 0) const = 0; virtual bool isFinite() const = 0; /// Returns `true` if `*this` is distributed or not. virtual bool isDistributed() const { return false; } protected: ID id; /// Underlying dof manager DOFManager & _dof_manager; UInt release_{0}; }; inline std::ostream & operator<<(std::ostream & stream, SolverVector & _this) { _this.printself(stream); return stream; } } // namespace akantu #endif /* AKANTU_SOLVER_VECTOR_HH_ */ diff --git a/src/solver/solver_vector_default.hh b/src/solver/solver_vector_default.hh index 5957ded5b..21cfa2d4f 100644 --- a/src/solver/solver_vector_default.hh +++ b/src/solver/solver_vector_default.hh @@ -1,157 +1,158 @@ /** * @file solver_vector_default.hh * * @author Nicolas Richart * * @date creation: Tue Jan 01 2019 * @date last modification: Sat May 23 2020 * * @brief Solver vector interface to Array * * * @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 "solver_vector.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SOLVER_VECTOR_DEFAULT_HH_ #define AKANTU_SOLVER_VECTOR_DEFAULT_HH_ namespace akantu { class DOFManagerDefault; } // namespace akantu namespace akantu { class SolverVectorArray : public SolverVector { public: SolverVectorArray(DOFManagerDefault & dof_manager, const ID & id); SolverVectorArray(const SolverVectorArray & vector, const ID & id); ~SolverVectorArray() override = default; virtual Array & getVector() = 0; virtual const Array & getVector() const = 0; void printself(std::ostream & stream, int indent = 0) const override { std::string space(indent, AKANTU_INDENT); stream << space << "SolverVectorArray [" << std::endl; stream << space << " + id: " << id << std::endl; this->getVector().printself(stream, indent + 1); stream << space << "]" << std::endl; } using SolverVector::isDistributed; }; /* -------------------------------------------------------------------------- */ template class SolverVectorArrayTmpl : public SolverVectorArray { public: SolverVectorArrayTmpl(DOFManagerDefault & dof_manager, Array_ & vector, const ID & id = "solver_vector_default") : SolverVectorArray(dof_manager, id), dof_manager(dof_manager), vector(vector) {} template ::value> * = nullptr> SolverVectorArrayTmpl(DOFManagerDefault & dof_manager, const ID & id = "solver_vector_default") : SolverVectorArray(dof_manager, id), dof_manager(dof_manager), vector(0, 1, id + ":vector") {} SolverVectorArrayTmpl(const SolverVectorArrayTmpl & vector, const ID & id = "solver_vector_default") : SolverVectorArray(vector, id), dof_manager(vector.dof_manager), vector(vector.vector) {} operator const Array &() const override { return getVector(); }; virtual operator Array &() { return getVector(); }; SolverVector & operator+=(const SolverVector & y) override; SolverVector & operator-=(const SolverVector & y) override; SolverVector & operator=(const SolverVector & y) override; void copy(const SolverVector & y) override; SolverVector & operator*=(const Real & alpha) override; void add(const SolverVector & y, const Real & alpha) override; Real dot(const SolverVector & y) const override; Real norm() const override; + Real normFreeDOFs() const override; void resize() override { static_assert(not std::is_const>::value, "Cannot resize a const Array"); this->vector.resize(this->localSize(), 0.); ++this->release_; } void set(Real val) override { static_assert(not std::is_const>::value, "Cannot clear a const Array"); this->vector.set(val); ++this->release_; } virtual bool isDistributed() const override { return false; } public: Array & getVector() override { return vector; } const Array & getVector() const override { return vector; } Int size() const override; Int localSize() const override; virtual Array & getGlobalVector() { return this->vector; } virtual void setGlobalVector(const Array & solution) { this->vector.copy(solution); } bool isFinite() const override { return vector.isFinite(); } protected: DOFManagerDefault & dof_manager; Array_ vector; template friend class SolverVectorArrayTmpl; }; /* -------------------------------------------------------------------------- */ using SolverVectorDefault = SolverVectorArrayTmpl>; /* -------------------------------------------------------------------------- */ template using SolverVectorDefaultWrap = SolverVectorArrayTmpl; template decltype(auto) make_solver_vector_default_wrap(DOFManagerDefault & dof_manager, Array & vector) { return SolverVectorDefaultWrap(dof_manager, vector); } } // namespace akantu /* -------------------------------------------------------------------------- */ #include "solver_vector_default_tmpl.hh" /* -------------------------------------------------------------------------- */ #endif /* AKANTU_SOLVER_VECTOR_DEFAULT_HH_ */ diff --git a/src/solver/solver_vector_default_tmpl.hh b/src/solver/solver_vector_default_tmpl.hh index 5857c7034..cc1453f4e 100644 --- a/src/solver/solver_vector_default_tmpl.hh +++ b/src/solver/solver_vector_default_tmpl.hh @@ -1,141 +1,156 @@ /** * @file solver_vector_default_tmpl.hh * * @author Nicolas Richart * * @date creation: Tue Jan 01 2019 * @date last modification: Sat May 23 2020 * * @brief Solver vector interface to Array * * * @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 "dof_manager_default.hh" #include "solver_vector_default.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SOLVER_VECTOR_DEFAULT_TMPL_HH_ #define AKANTU_SOLVER_VECTOR_DEFAULT_TMPL_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ inline SolverVectorArray::SolverVectorArray(DOFManagerDefault & dof_manager, const ID & id) : SolverVector(dof_manager, id) {} /* -------------------------------------------------------------------------- */ inline SolverVectorArray::SolverVectorArray(const SolverVectorArray & vector, const ID & id) : SolverVector(vector, id) {} /* -------------------------------------------------------------------------- */ template SolverVector & SolverVectorArrayTmpl::operator+=(const SolverVector & y) { const auto & y_ = aka::as_type(y); this->vector += y_.getVector(); ++this->release_; return *this; } /* -------------------------------------------------------------------------- */ template SolverVector & SolverVectorArrayTmpl::operator-=(const SolverVector & y) { const auto & y_ = aka::as_type(y); this->vector -= y_.getVector(); ++this->release_; return *this; } /* -------------------------------------------------------------------------- */ template SolverVector & SolverVectorArrayTmpl::operator=(const SolverVector & y) { const auto & y_ = aka::as_type(y); this->vector.copy(y_.getVector()); this->release_ = y.release(); return *this; } /* -------------------------------------------------------------------------- */ template void SolverVectorArrayTmpl::copy(const SolverVector & y) { this->operator=(y); } /* -------------------------------------------------------------------------- */ template SolverVector & SolverVectorArrayTmpl::operator*=(const Real & alpha) { this->vector *= alpha; ++this->release_; return *this; } /* -------------------------------------------------------------------------- */ template void SolverVectorArrayTmpl::add(const SolverVector & y, const Real & alpha) { const auto & y_ = aka::as_type(y); for (auto && data : zip(this->vector, y_.getVector())) { auto && x = std::get<0>(data); auto && y = std::get<1>(data); x += y * alpha; } ++this->release_; } /* -------------------------------------------------------------------------- */ template Real SolverVectorArrayTmpl::dot(const SolverVector & y) const { const auto & y_ = aka::as_type(y); Real sum = 0.; for (auto [d, a, b] : enumerate(this->vector, y_.getVector())) { sum += this->dof_manager.isLocalOrMasterDOF(d) * a * b; } return sum; } /* -------------------------------------------------------------------------- */ template Real SolverVectorArrayTmpl::norm() const { auto a = this->dot(*this); return std::sqrt(a); } +/* -------------------------------------------------------------------------- */ +template +Real SolverVectorArrayTmpl::normFreeDOFs() const { + const auto & blocked_dofs = this->dof_manager.getBlockedDOFs(); + + Real sum = 0.; + for (auto [d, c, a] : enumerate(blocked_dofs, this->vector)) { + bool is_local_node = this->dof_manager.isLocalOrMasterDOF(d); + bool is_blocked_dof = c; + if ((!is_blocked_dof) && is_local_node) { + sum += a * a; + } + } + return std::sqrt(sum); +} /* -------------------------------------------------------------------------- */ template inline Int SolverVectorArrayTmpl::size() const { return this->dof_manager.getSystemSize(); } /* -------------------------------------------------------------------------- */ template inline Int SolverVectorArrayTmpl::localSize() const { return dof_manager.getLocalSystemSize(); } } // namespace akantu #endif /* AKANTU_SOLVER_VECTOR_DEFAULT_TMPL_HH_ */ diff --git a/src/solver/solver_vector_distributed.cc b/src/solver/solver_vector_distributed.cc index 926a3f911..d8c31ffb4 100644 --- a/src/solver/solver_vector_distributed.cc +++ b/src/solver/solver_vector_distributed.cc @@ -1,99 +1,107 @@ /** * @file solver_vector_distributed.cc * * @author Nicolas Richart * * @date creation: Tue Jan 01 2019 * @date last modification: Sat May 23 2020 * * @brief Solver vector interface for distributed arrays * * * @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 "solver_vector_distributed.hh" #include "dof_manager_default.hh" #include "dof_synchronizer.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ SolverVectorDistributed::SolverVectorDistributed( DOFManagerDefault & dof_manager, const ID & id) : SolverVectorDefault(dof_manager, id) {} /* -------------------------------------------------------------------------- */ SolverVectorDistributed::SolverVectorDistributed( const SolverVectorDefault & vector, const ID & id) : SolverVectorDefault(vector, id) {} /* -------------------------------------------------------------------------- */ Array & SolverVectorDistributed::getGlobalVector() { auto & synchronizer = dof_manager.getSynchronizer(); if (not this->global_vector) { this->global_vector = std::make_unique>(0, 1, "global_residual"); } if (synchronizer.getCommunicator().whoAmI() == 0) { this->global_vector->resize(dof_manager.getSystemSize()); synchronizer.gather(this->vector, *this->global_vector); } else { synchronizer.gather(this->vector); } return *this->global_vector; } /* -------------------------------------------------------------------------- */ Real SolverVectorDistributed::dot(const SolverVector & y) const { auto sum = this->SolverVectorDefault::dot(y); auto & synchronizer = dof_manager.getSynchronizer(); synchronizer.getCommunicator().allReduce(sum, SynchronizerOperation::_sum); return sum; } +/* -------------------------------------------------------------------------- */ +Real SolverVectorDistributed::normFreeDOFs() const { + + auto sum = this->SolverVectorDefault::normFreeDOFs(); + auto & synchronizer = dof_manager.getSynchronizer(); + synchronizer.getCommunicator().allReduce(sum, SynchronizerOperation::_sum); + return sum; +} /* -------------------------------------------------------------------------- */ void SolverVectorDistributed::setGlobalVector(const Array & solution) { auto & synchronizer = dof_manager.getSynchronizer(); if (synchronizer.getCommunicator().whoAmI() == 0) { synchronizer.scatter(this->vector, solution); } else { synchronizer.scatter(this->vector); } } /* -------------------------------------------------------------------------- */ bool SolverVectorDistributed::isFinite() const { auto & synchronizer = dof_manager.getSynchronizer(); bool is_finite = this->vector.isFinite(); synchronizer.getCommunicator().allReduce(is_finite, SynchronizerOperation::_land); return is_finite; } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/solver/solver_vector_distributed.hh b/src/solver/solver_vector_distributed.hh index d1aa385ba..a2dde59f6 100644 --- a/src/solver/solver_vector_distributed.hh +++ b/src/solver/solver_vector_distributed.hh @@ -1,65 +1,67 @@ /** * @file solver_vector_distributed.hh * * @author Nicolas Richart * * @date creation: Thu Feb 21 2013 * @date last modification: Tue Jan 01 2019 * * @brief Solver vector interface for distributed arrays * * * @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 "solver_vector_default.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_SOLVER_VECTOR_DISTRIBUTED_HH_ #define AKANTU_SOLVER_VECTOR_DISTRIBUTED_HH_ namespace akantu { class SolverVectorDistributed : public SolverVectorDefault { public: SolverVectorDistributed(DOFManagerDefault & dof_manager, const ID & id = "solver_vector_mumps"); SolverVectorDistributed(const SolverVectorDefault & vector, const ID & id = "solver_vector_mumps"); Array & getGlobalVector() override; void setGlobalVector(const Array & solution) override; bool isFinite() const override; Real dot(const SolverVector & y) const override; + Real normFreeDOFs() const override; + virtual bool isDistributed() const override { return true; } protected: // full vector in case it needs to be centralized on master std::unique_ptr> global_vector; }; } // namespace akantu #endif /* AKANTU_SOLVER_VECTOR_DISTRIBUTED_HH_ */