diff --git a/python/py_aka_common.cc b/python/py_aka_common.cc index 7a5ef2c77..07fbcc881 100644 --- a/python/py_aka_common.cc +++ b/python/py_aka_common.cc @@ -1,111 +1,112 @@ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; namespace akantu { /* -------------------------------------------------------------------------- */ #define PY_AKANTU_PP_VALUE(s, data, elem) \ .value(BOOST_PP_STRINGIZE(elem), BOOST_PP_CAT(data, elem)) #define PY_AKANTU_REGISTER_ENUM_(type_name, list, prefix, mod) \ py::enum_(mod, BOOST_PP_STRINGIZE(type_name)) \ BOOST_PP_SEQ_FOR_EACH(PY_AKANTU_PP_VALUE, prefix, list) \ .export_values() #define PY_AKANTU_REGISTER_CLASS_ENUM(type_name, list, mod) \ PY_AKANTU_REGISTER_ENUM_(type_name, list, type_name::_, mod) #define PY_AKANTU_REGISTER_ENUM(type_name, list, mod) \ PY_AKANTU_REGISTER_ENUM_(type_name, list, , mod) /* -------------------------------------------------------------------------- */ void register_initialize(py::module & mod) { mod.def("__initialize", []() { int nb_args = 0; char ** null = nullptr; initialize(nb_args, null); }); } void register_enums(py::module & mod) { py::enum_(mod, "SpatialDirection") .value("_x", _x) .value("_y", _y) .value("_z", _z) .export_values(); py::enum_(mod, "AnalysisMethod") .value("_static", _static) .value("_implicit_dynamic", _implicit_dynamic) .value("_explicit_lumped_mass", _explicit_lumped_mass) .value("_explicit_lumped_capacity", _explicit_lumped_capacity) .value("_explicit_consistent_mass", _explicit_consistent_mass) .export_values(); + PY_AKANTU_REGISTER_CLASS_ENUM(ModelType, AKANTU_MODEL_TYPES, mod); PY_AKANTU_REGISTER_CLASS_ENUM(NonLinearSolverType, AKANTU_NON_LINEAR_SOLVER_TYPES, mod); PY_AKANTU_REGISTER_CLASS_ENUM(TimeStepSolverType, AKANTU_TIME_STEP_SOLVER_TYPE, mod); PY_AKANTU_REGISTER_CLASS_ENUM(IntegrationSchemeType, AKANTU_INTEGRATION_SCHEME_TYPE, mod); PY_AKANTU_REGISTER_CLASS_ENUM(SolveConvergenceCriteria, AKANTU_SOLVE_CONVERGENCE_CRITERIA, mod); py::enum_(mod, "CohesiveMethod") .value("_intrinsic", _intrinsic) .value("_extrinsic", _extrinsic) .export_values(); py::enum_(mod, "GhostType") .value("_not_ghost", _not_ghost) .value("_ghost", _ghost) .value("_casper", _casper) .export_values(); py::enum_(mod, "MeshIOType") .value("_miot_auto", _miot_auto) .value("_miot_gmsh", _miot_gmsh) .value("_miot_gmsh_struct", _miot_gmsh_struct) .value("_miot_diana", _miot_diana) .value("_miot_abaqus", _miot_abaqus) .export_values(); py::enum_(mod, "MatrixType") .value("_unsymmetric", _unsymmetric) .value("_symmetric", _symmetric) .export_values(); PY_AKANTU_REGISTER_ENUM(ElementType, AKANTU_ALL_ELEMENT_TYPE(_not_defined), mod); PY_AKANTU_REGISTER_ENUM(ElementKind, AKANTU_ELEMENT_KIND(_ek_not_defined), mod); } /* -------------------------------------------------------------------------- */ #define AKANTU_PP_STR_TO_TYPE2(s, data, elem) ({BOOST_PP_STRINGIZE(elem), elem}) void register_functions(py::module & mod) { mod.def("getElementTypes", []() { std::map element_types{ BOOST_PP_SEQ_FOR_EACH_I( AKANTU_PP_ENUM, BOOST_PP_SEQ_SIZE(AKANTU_ek_regular_ELEMENT_TYPE), BOOST_PP_SEQ_TRANSFORM(AKANTU_PP_STR_TO_TYPE2, akantu, AKANTU_ek_regular_ELEMENT_TYPE))}; return element_types; }); } #undef AKANTU_PP_STR_TO_TYPE2 } // namespace akantu diff --git a/python/py_model.cc b/python/py_model.cc index 7f1c35e28..7d9a63e34 100644 --- a/python/py_model.cc +++ b/python/py_model.cc @@ -1,81 +1,89 @@ /* -------------------------------------------------------------------------- */ #include "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ void register_model(py::module & mod) { py::class_(mod, "DOFManager") .def("getMatrix", &DOFManager::getMatrix, py::return_value_policy::reference) .def( "getNewMatrix", [](DOFManager & self, const std::string & name, const std::string & matrix_to_copy_id) -> decltype(auto) { return self.getNewMatrix(name, matrix_to_copy_id); }, py::return_value_policy::reference) .def( "getResidual", [](DOFManager & self) -> decltype(auto) { return self.getResidual(); }, py::return_value_policy::reference) .def("getArrayPerDOFs", &DOFManager::getArrayPerDOFs) .def("assembleToResidual", &DOFManager::assembleToResidual); py::class_(mod, "NonLinearSolver") .def( "set", [](NonLinearSolver & self, const std::string & id, const Real & val) { if (id == "max_iterations") { self.set(id, int(val)); } else { self.set(id, val); } }) .def("set", [](NonLinearSolver & self, const std::string & id, const SolveConvergenceCriteria & val) { self.set(id, val); }); py::class_(mod, "ModelSolver", py::multiple_inheritance()) .def("getNonLinearSolver", (NonLinearSolver & (ModelSolver::*)(const ID &)) & ModelSolver::getNonLinearSolver, py::arg("solver_id") = "", py::return_value_policy::reference) .def("solveStep", [](ModelSolver & self) { self.solveStep(); }) .def("solveStep", [](ModelSolver & self, const ID & solver_id) { self.solveStep(solver_id); }); py::class_(mod, "Model", py::multiple_inheritance()) .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("addDumpFieldVectorToDumper", &Model::addDumpFieldVectorToDumper) .def("addDumpFieldToDumper", &Model::addDumpFieldToDumper) .def("dump", &Model::dump) .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) { + self.setIntegrationScheme(id, primal, scheme); + }) .def("getDOFManager", &Model::getDOFManager, py::return_value_policy::reference); } } // namespace akantu diff --git a/python/py_phase_field_model.cc b/python/py_phase_field_model.cc index b06c9609a..ac7866866 100644 --- a/python/py_phase_field_model.cc +++ b/python/py_phase_field_model.cc @@ -1,122 +1,124 @@ /* -------------------------------------------------------------------------- */ #include "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #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, \ [](PhaseFieldModel & self) -> decltype(auto) { \ return self.func_name(); \ }, \ py::return_value_policy::reference) #define def_function(func_name) \ def(#func_name, [](PhaseFieldModel & self) -> decltype(auto) { \ return self.func_name(); \ }) /* -------------------------------------------------------------------------- */ [[gnu::visibility("default")]] void register_phase_field_model(py::module & mod) { py::class_(mod, "PhaseFieldModelOptions") .def(py::init(), py::arg("analysis_method") = _static); py::class_(mod, "PhaseFieldModel", py::multiple_inheritance()) .def(py::init(), py::arg("mesh"), py::arg("spatial_dimension") = _all_dimensions, py::arg("id") = "phase_field_model", py::arg("memory_id") = 0, py::arg("model_type") = ModelType::_phase_field_model) .def("initFull", [](PhaseFieldModel & self, const PhaseFieldModelOptions & options) { self.initFull(options); }, py::arg("_analysis_method") = PhaseFieldModelOptions()) .def("initFull", [](PhaseFieldModel & self, const AnalysisMethod & analysis_method) { self.initFull(_analysis_method = analysis_method); }, py::arg("_analysis_method")) .def_deprecated("applyDirichletBC", "Deprecated: use applyBC") .def("applyBC", [](PhaseFieldModel & self, BC::Dirichlet::DirichletFunctor & func, const std::string & element_group) { self.applyBC(func, element_group); }) .def("applyBC", [](PhaseFieldModel & self, BC::Neumann::NeumannFunctor & func, const std::string & element_group) { self.applyBC(func, element_group); }) .def("setTimeStep", &PhaseFieldModel::setTimeStep, py::arg("time_step"), py::arg("solver_id") = "") .def_function(assembleStiffnessMatrix) .def_function(assembleInternalForces) .def_function_nocopy(getDamage) .def_function_nocopy(getInternalForce) .def_function_nocopy(getBlockedDOFs) .def_function_nocopy(getMesh) .def("dump", py::overload_cast<>(&PhaseFieldModel::dump)) .def("dump", py::overload_cast(&PhaseFieldModel::dump)) .def("dump", py::overload_cast( &PhaseFieldModel::dump)) .def("dump", py::overload_cast( &PhaseFieldModel::dump)) .def("getPhaseField", py::overload_cast(&PhaseFieldModel::getPhaseField), py::return_value_policy::reference) .def("getPhaseField", py::overload_cast( &PhaseFieldModel::getPhaseField), py::return_value_policy::reference) .def("getPhaseFieldIndex", &PhaseFieldModel::getPhaseFieldIndex) .def("setPhaseFieldSelector", &PhaseFieldModel::setPhaseFieldSelector); } [[gnu::visibility("default")]] void register_phase_field_coupler(py::module & mod) { py::class_(mod, "CouplerSolidPhaseField") .def(py::init(), py::arg("mesh"), py::arg("spatial_dimension") = _all_dimensions, py::arg("id") = "coupler_solid_phasefield", py::arg("memory_id") = 0, py::arg("model_type") = ModelType::_coupler_solid_phasefield) - .def("solve", [](CouplerSolidPhaseField & self, const ID & solver_id) { - self.solve(solver_id); }) + .def("solve", [](CouplerSolidPhaseField & self, const ID & solid_solver_id, + const ID & phase_solver_id) { + self.solve(solid_solver_id, phase_solver_id); + }) .def("getSolidMechanicsModel", &CouplerSolidPhaseField::getSolidMechanicsModel, py::return_value_policy::reference) .def("getPhaseFieldModel", &CouplerSolidPhaseField::getPhaseFieldModel, py::return_value_policy::reference); } } diff --git a/src/model/model_couplers/coupler_solid_phasefield.cc b/src/model/model_couplers/coupler_solid_phasefield.cc index 8c6761b74..db74b7f10 100644 --- a/src/model/model_couplers/coupler_solid_phasefield.cc +++ b/src/model/model_couplers/coupler_solid_phasefield.cc @@ -1,704 +1,709 @@ /** * @file coupler_solid_phasefield.cc * * @author Mohit Pundir * * @date creation: Fri Sep 28 2018 * @date last modification: Thu Jun 20 2019 * * @brief class for coupling of solid mechancis and phase model * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "coupler_solid_phasefield.hh" #include "dumpable_inline_impl.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #include "element_synchronizer.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { CouplerSolidPhaseField::CouplerSolidPhaseField(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id, const ModelType model_type) : Model(mesh, model_type, dim, id, memory_id) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("CouplerSolidPhaseField", mesh, Model::spatial_dimension); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("coupler_solid_phasefield", id, true); this->mesh.addDumpMeshToDumper("coupler_solid_phasefield", mesh, Model::spatial_dimension, _not_ghost, _ek_regular); #endif this->registerDataAccessor(*this); solid = new SolidMechanicsModel(mesh, Model::spatial_dimension, "solid_mechanics_model", 0); phase = new PhaseFieldModel(mesh, Model::spatial_dimension, "phase_field_model", 0); if (this->mesh.isDistributed()) { auto & synchronizer = this->mesh.getElementSynchronizer(); this->registerSynchronizer(synchronizer, SynchronizationTag::_csp_damage); this->registerSynchronizer(synchronizer, SynchronizationTag::_csp_strain); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ CouplerSolidPhaseField::~CouplerSolidPhaseField() {} /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); this->initBC(*this, *displacement, *displacement_increment, *external_force); solid->initFull( _analysis_method = this->method); phase->initFull( _analysis_method = this->method); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::initModel() { getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ FEEngine & CouplerSolidPhaseField::getFEEngineBoundary(const ID & name) { return dynamic_cast( getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) { auto & solid_model_solver = aka::as_type(*solid); solid_model_solver.initSolver(time_step_solver_type, non_linear_solver_type); auto & phase_model_solver = aka::as_type(*phase); phase_model_solver.initSolver(time_step_solver_type, non_linear_solver_type); } /* -------------------------------------------------------------------------- */ std::tuple CouplerSolidPhaseField::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_lumped_mass: { return std::make_tuple("explicit_lumped", TimeStepSolverType::_dynamic_lumped); } case _explicit_consistent_mass: { return std::make_tuple("explicit", TimeStepSolverType::_dynamic); } case _static: { return std::make_tuple("static", TimeStepSolverType::_static); } case _implicit_dynamic: { return std::make_tuple("implicit", TimeStepSolverType::_dynamic); } default: return std::make_tuple("unknown", TimeStepSolverType::_not_defined); } } /* -------------------------------------------------------------------------- */ TimeStepSolverType CouplerSolidPhaseField::getDefaultSolverType() const { return TimeStepSolverType::_dynamic_lumped; } /* -------------------------------------------------------------------------- */ ModelSolverOptions CouplerSolidPhaseField::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_dynamic_lumped: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_dynamic: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_static: { options.non_linear_solver_type = NonLinearSolverType::_linear; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); break; } return options; } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleResidual() { // computes the internal forces this->assembleInternalForces(); auto & solid_internal_force = solid->getInternalForce(); auto & solid_external_force = solid->getExternalForce(); auto & phasefield_internal_force = phase->getInternalForce(); auto & phasefield_external_force = phase->getExternalForce(); /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", solid_external_force, 1); this->getDOFManager().assembleToResidual("displacement", solid_internal_force, 1); this->getDOFManager().assembleToResidual("damage", phasefield_external_force, 1); this->getDOFManager().assembleToResidual("damage", phasefield_internal_force, 1); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleResidual(const ID & residual_part) { AKANTU_DEBUG_IN(); auto & solid_internal_force = solid->getInternalForce(); auto & solid_external_force = solid->getExternalForce(); auto & phasefield_internal_force = phase->getInternalForce(); auto & phasefield_external_force = phase->getExternalForce(); if ("external" == residual_part) { this->getDOFManager().assembleToResidual("displacement", solid_external_force, 1); this->getDOFManager().assembleToResidual("displacement", solid_internal_force, 1); AKANTU_DEBUG_OUT(); return; } if ("internal" == residual_part) { this->getDOFManager().assembleToResidual("damage", phasefield_external_force, 1); this->getDOFManager().assembleToResidual("damage", phasefield_internal_force, 1); AKANTU_DEBUG_OUT(); return; } AKANTU_CUSTOM_EXCEPTION( debug::SolverCallbackResidualPartUnknown(residual_part)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::predictor() { auto & solid_model_solver = aka::as_type(*solid); solid_model_solver.predictor(); auto & phase_model_solver = aka::as_type(*phase); phase_model_solver.predictor(); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::corrector() { auto & solid_model_solver = aka::as_type(*solid); solid_model_solver.corrector(); auto & phase_model_solver = aka::as_type(*phase); phase_model_solver.corrector(); } /* -------------------------------------------------------------------------- */ MatrixType CouplerSolidPhaseField::getMatrixType(const ID & matrix_id) { if (matrix_id == "K") return _symmetric; if (matrix_id == "M") { return _symmetric; } return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } else if (matrix_id == "M") { solid->assembleMass(); } } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleLumpedMatrix(const ID & matrix_id) { if (matrix_id == "M") { solid->assembleMassLumped(); } } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::beforeSolveStep() { auto & solid_solver_callback = aka::as_type(*solid); solid_solver_callback.beforeSolveStep(); auto & phase_solver_callback = aka::as_type(*phase); phase_solver_callback.beforeSolveStep(); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::afterSolveStep(bool converged) { auto & solid_solver_callback = aka::as_type(*solid); solid_solver_callback.afterSolveStep(converged); auto & phase_solver_callback = aka::as_type(*phase); phase_solver_callback.afterSolveStep(converged); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); solid->assembleInternalForces(); phase->assembleInternalForces(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix"); solid->assembleStiffnessMatrix(); phase->assembleStiffnessMatrix(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleMassLumped() { solid->assembleMassLumped(); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleMass() { solid->assembleMass(); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleMassLumped(GhostType ghost_type) { solid->assembleMassLumped(ghost_type); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::assembleMass(GhostType ghost_type) { solid->assembleMass(ghost_type); } /* ------------------------------------------------------------------------- */ void CouplerSolidPhaseField::computeDamageOnQuadPoints( const GhostType & ghost_type) { AKANTU_DEBUG_IN(); auto & fem = phase->getFEEngine(); auto & mesh = phase->getMesh(); auto nb_materials = solid->getNbMaterials(); auto nb_phasefields = phase->getNbPhaseFields(); AKANTU_DEBUG_ASSERT(nb_phasefields == nb_materials, "The number of phasefields and materials should be equal" ); for(auto index : arange(nb_materials)) { auto & material = solid->getMaterial(index); for(auto index2 : arange(nb_phasefields)) { auto & phasefield = phase->getPhaseField(index2); if(phasefield.getName() == material.getName()){ switch (spatial_dimension) { case 1: { auto & mat = static_cast &>(material); auto & damage = mat.getDamage(); for (auto & type : mesh.elementTypes(Model::spatial_dimension, ghost_type)) { auto & damage_on_qpoints_vect = damage(type, ghost_type); fem.interpolateOnIntegrationPoints(phase->getDamage(), damage_on_qpoints_vect, 1, type, ghost_type); } - break; } case 2: { auto & mat = static_cast &>(material); auto & damage = mat.getDamage(); for (auto & type : mesh.elementTypes(Model::spatial_dimension, ghost_type)) { auto & damage_on_qpoints_vect = damage(type, ghost_type); fem.interpolateOnIntegrationPoints(phase->getDamage(), damage_on_qpoints_vect, 1, type, ghost_type); } break; } default: auto & mat = static_cast &>(material); + auto & damage = mat.getDamage(); + + for (auto & type : + mesh.elementTypes(Model::spatial_dimension, ghost_type)) { + auto & damage_on_qpoints_vect = damage(type, ghost_type); + fem.interpolateOnIntegrationPoints(phase->getDamage(), damage_on_qpoints_vect, + 1, type, ghost_type); + } break; } } } } AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------- */ void CouplerSolidPhaseField::computeStrainOnQuadPoints( const GhostType & ghost_type) { AKANTU_DEBUG_IN(); auto & mesh = solid->getMesh(); auto nb_materials = solid->getNbMaterials(); auto nb_phasefields = phase->getNbPhaseFields(); AKANTU_DEBUG_ASSERT(nb_phasefields == nb_materials, "The number of phasefields and materials should be equal" ); for(auto index : arange(nb_materials)) { auto & material = solid->getMaterial(index); for(auto index2 : arange(nb_phasefields)) { auto & phasefield = phase->getPhaseField(index2); if(phasefield.getName() == material.getName()){ auto & strain_on_qpoints = phasefield.getStrain(); auto & gradu_on_qpoints = material.getGradU(); for (auto & type: mesh.elementTypes(spatial_dimension, ghost_type)) { auto & strain_on_qpoints_vect = strain_on_qpoints(type, ghost_type); auto & gradu_on_qpoints_vect = gradu_on_qpoints(type, ghost_type); for (auto && values: zip(make_view(strain_on_qpoints_vect, spatial_dimension, spatial_dimension), make_view(gradu_on_qpoints_vect, spatial_dimension, spatial_dimension))) { auto & strain = std::get<0>(values); auto & grad_u = std::get<1>(values); gradUToEpsilon(grad_u, strain); } } break; } } } AKANTU_DEBUG_OUT(); } /* ------------------------------------------------------------------------- */ -void CouplerSolidPhaseField::solve(const ID & solver_id) { +void CouplerSolidPhaseField::solve(const ID & solid_solver_id, const ID & phase_solver_id) { - solid->solveStep(solver_id); + solid->solveStep(solid_solver_id); this->computeStrainOnQuadPoints(_not_ghost); - phase->solveStep(); + phase->solveStep(phase_solver_id); this->computeDamageOnQuadPoints(_not_ghost); solid->assembleInternalForces(); } /* ------------------------------------------------------------------------- */ void CouplerSolidPhaseField::gradUToEpsilon(const Matrix & grad_u, Matrix & epsilon) { for (UInt i = 0; i < Model::spatial_dimension; ++i) { for (UInt j = 0; j < Model::spatial_dimension; ++j) epsilon(i, j) = 0.5 * (grad_u(i, j) + grad_u(j, i)); } } /* ------------------------------------------------------------------------- */ bool CouplerSolidPhaseField::checkConvergence(Array & u_new, Array & u_old, Array & d_new, Array & d_old) { const Array & blocked_dofs = solid->getBlockedDOFs(); UInt nb_degree_of_freedom = u_new.size(); auto u_n_it = u_new.begin(); auto u_o_it = u_old.begin(); auto bld_it = blocked_dofs.begin(); Real norm = 0; for (UInt n = 0; n < nb_degree_of_freedom; ++n, ++u_n_it, ++u_o_it, ++bld_it) { if ((!*bld_it)) { norm += (*u_n_it - *u_o_it) * (*u_n_it - *u_o_it); } } norm = std::sqrt(norm); auto d_n_it = d_new.begin(); auto d_o_it = d_old.begin(); nb_degree_of_freedom = d_new.size(); Real norm2 = 0; for (UInt i = 0; i < nb_degree_of_freedom; ++i) { norm2 += (*d_n_it - *d_o_it); } norm2 = std::sqrt(norm2); Real error = std::max(norm, norm2); Real tolerance = 1e-8; if (error < tolerance) { return true; } return false; } /* -------------------------------------------------------------------------- */ /*UInt CouplerSolidPhaseField::getNbData(const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; switch (tag) { case SynchronizationTag::_csp_damage: { size += getNbIntegrationPoints(elements) * sizeof(Real); break; } case SynchronizationTag::_csp_strain: { size += getNbIntegrationPoints(elements) * spatial_dimension * spatial_dimension * sizeof(Real); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } -/* -------------------------------------------------------------------------- */ -/*void CouplerSolidPhaseField::packData(CommunicationBuffer & buffer, +void CouplerSolidPhaseField::packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { auto & mat = static_cast &>(solid->getMaterial(0)); switch (tag) { case SynchronizationTag::_csp_damage: { packElementalDataHelper(mat.getDamage(), buffer, elements, true, getFEEngine()); break; } case SynchronizationTag::_csp_strain: { packElementalDataHelper(phase->getStrain(), buffer, elements, true, getFEEngine()); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } - }*/ + } -/* -------------------------------------------------------------------------- */ -/*void CouplerSolidPhaseField::unpackData(CommunicationBuffer & buffer, +void CouplerSolidPhaseField::unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { auto & mat = static_cast &>(solid->getMaterial(0)); switch (tag) { case SynchronizationTag::_csp_damage: { unpackElementalDataHelper(mat.getDamage(), buffer, elements, true, getFEEngine()); break; } case SynchronizationTag::_csp_strain: { unpackElementalDataHelper(phase->getStrain(), buffer, elements, true, getFEEngine()); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } }*/ /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidPhaseField::createElementalField( const std::string & field_name, const std::string & group_name, bool padding_flag, UInt spatial_dimension, ElementKind kind) { return solid->createElementalField(field_name, group_name, padding_flag, spatial_dimension, kind); std::shared_ptr field; return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidPhaseField::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { return solid->createNodalFieldReal(field_name, group_name, padding_flag); std::shared_ptr field; return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidPhaseField::createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) { return solid->createNodalFieldBool(field_name, group_name, padding_flag); std::shared_ptr field; return field; } #else /* -------------------------------------------------------------------------- */ std::shared_ptr CouplerSolidPhaseField::createElementalField( const std::string &, const std::string &, bool, UInt , ElementKind) { return nullptr; } /* ----------------------------------------------------------------------- */ std::shared_ptr CouplerSolidPhaseField::createNodalFieldReal(const std::string &, const std::string &, bool) { return nullptr; } /*-------------------------------------------------------------------*/ std::shared_ptr CouplerSolidPhaseField::createNodalFieldBool(const std::string &, const std::string &, bool) { return nullptr; } #endif /* -----------------------------------------------------------------------*/ void CouplerSolidPhaseField::dump(const std::string & dumper_name) { solid->onDump(); mesh.dump(dumper_name); } /* ------------------------------------------------------------------------*/ void CouplerSolidPhaseField::dump(const std::string & dumper_name, UInt step) { solid->onDump(); mesh.dump(dumper_name, step); } /* ----------------------------------------------------------------------- */ void CouplerSolidPhaseField::dump(const std::string & dumper_name, Real time, UInt step) { solid->onDump(); mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::dump() { solid->onDump(); mesh.dump(); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::dump(UInt step) { solid->onDump(); mesh.dump(step); } /* -------------------------------------------------------------------------- */ void CouplerSolidPhaseField::dump(Real time, UInt step) { solid->onDump(); mesh.dump(time, step); } } // namespace akantu diff --git a/src/model/model_couplers/coupler_solid_phasefield.hh b/src/model/model_couplers/coupler_solid_phasefield.hh index 01f3aae09..14a38c1ca 100644 --- a/src/model/model_couplers/coupler_solid_phasefield.hh +++ b/src/model/model_couplers/coupler_solid_phasefield.hh @@ -1,289 +1,292 @@ /** * @file solid_phase_coupler.hh * * @author Mohit Pundir * * @date creation: Fri Sep 28 2018 * @date last modification: Fri Sep 28 2018 * * @brief class for coupling of solid mechancis and phasefield model * * @section LICENSE * * Copyright (©) 2010-2018 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 "material.hh" #include "material_phasefield.hh" #include "model.hh" #include "phase_field_model.hh" #include "solid_mechanics_model.hh" #include "sparse_matrix.hh" #include "time_step_solver.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COUPLER_SOLID_PHASEFIELD_HH__ #define __AKANTU_COUPLER_SOLID_PHASEFIELD_HH__ /* ------------------------------------------------------------------------ */ /* Coupling : Solid Mechanics / PhaseField */ /* ------------------------------------------------------------------------ */ namespace akantu { template class IntegratorGauss; template class ShapeLagrange; class DOFManager; } // namespace akantu namespace akantu { class CouplerSolidPhaseField : public Model, public DataAccessor, public DataAccessor, public BoundaryCondition { /* ------------------------------------------------------------------------ */ /* Constructor/Destructors */ /* ------------------------------------------------------------------------ */ using MyFEEngineType = FEEngineTemplate; public: CouplerSolidPhaseField( Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "coupler_solid_phasefield", const MemoryID & memory_id = 0, const ModelType model_type = ModelType::_coupler_solid_phasefield); ~CouplerSolidPhaseField() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize the complete model void initFullImpl(const ModelOptions & options) override; /// initialize the modelType void initModel() override; /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; /* ------------------------------------------------------------------------ */ /* Solver Interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the contact stiffness matrix virtual void assembleStiffnessMatrix(); /// assembles the contant internal forces virtual void assembleInternalForces(); public: /// computes damage on quad points for solid mechanics model from /// damage array from phasefield model void computeDamageOnQuadPoints(const GhostType &); /// computes strain on quadrature points for phasefield model from /// displacement gradient from solid mechanics model void computeStrainOnQuadPoints(const GhostType & ghost_type); /// solve the coupled model - void solve(const ID & solver_id = ""); + void solve(const ID & solid_solver_id = "", const ID & phase_solver_id = ""); private: /// computes small strain from displacement gradient void gradUToEpsilon(const Matrix & grad_u, Matrix & epsilon); /// test the convergence criteria bool checkConvergence(Array &, Array &, Array &, Array &); 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() override { return true; } /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) override; /// callback for the solver, this assembles different matrices void assembleMatrix(const ID & matrix_id) override; /// callback for the solver, this assembles the stiffness matrix void assembleLumpedMatrix(const ID & matrix_id) override; /// callback for the model to instantiate the matricess when needed void initSolver(TimeStepSolverType, NonLinearSolverType) override; /// callback for the solver, this is called at beginning of solve void predictor() override; /// callback for the solver, this is called at end of solve void corrector() override; /// callback for the solver, this is called at beginning of solve void beforeSolveStep() override; /// callback for the solver, this is called at end of solve void afterSolveStep(bool converged = true) override; /// solve the coupled model //void solveStep(const ID & solver_id = "") override; /// solve a step using a given pre instantiated time step solver and /// non linear solver with a user defined callback instead of the /// model itself /!\ This can mess up everything //void solveStep(SolverCallback & callback, const ID & solver_id = "") override; /* ------------------------------------------------------------------------ */ /* Mass matrix for solid mechanics model */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the mass matrix for consistent mass resolutions void assembleMass(); protected: /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumped(GhostType ghost_type); /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMass(GhostType ghost_type); protected: /* -------------------------------------------------------------------------- */ TimeStepSolverType getDefaultSolverType() const override; /* -------------------------------------------------------------------------- */ ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const; public: bool isDefaultSolverExplicit() { return method == _explicit_lumped_mass; } /* ------------------------------------------------------------------------ */ public: // DataAccessor + UInt getNbData(const Array &, const SynchronizationTag &) const override { return 0; } void packData(CommunicationBuffer &, const Array &, const SynchronizationTag &) const override {} void unpackData(CommunicationBuffer &, const Array &, - const SynchronizationTag &) override {} - - UInt getNbData(const Array & indexes, - const SynchronizationTag & tag) const override { + const SynchronizationTag &) override {} + + UInt getNbData(__attribute__((unused)) const Array & indexes, + __attribute__((unused)) const SynchronizationTag & tag) const override { return 0; } - void packData(CommunicationBuffer & buffer, const Array & dofs, - const SynchronizationTag & tag) const override{} - - void unpackData(CommunicationBuffer & buffer, const Array & dofs, - const SynchronizationTag & tag) override {} + void packData(__attribute__((unused)) CommunicationBuffer & buffer, + __attribute__((unused)) const Array & dofs, + __attribute__((unused)) const SynchronizationTag & tag) const override{} + void unpackData(__attribute__((unused)) CommunicationBuffer & buffer, + __attribute__((unused)) const Array & dofs, + __attribute__((unused)) const SynchronizationTag & tag) override {} + /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: FEEngine & getFEEngineBoundary(const ID & name = "") override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt); /// get the solid mechanics model AKANTU_GET_MACRO(SolidMechanicsModel, *solid, SolidMechanicsModel &); /// get the contact mechanics model AKANTU_GET_MACRO(PhaseFieldModel, *phase, PhaseFieldModel &); /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: std::shared_ptr createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, UInt spatial_dimension, ElementKind kind) override; virtual void dump(const std::string & dumper_name); virtual void dump(const std::string & dumper_name, UInt step); virtual void dump(const std::string & dumper_name, Real time, UInt step); void dump() override; virtual void dump(UInt step); virtual void dump(Real time, UInt step); /* ------------------------------------------------------------------------ */ /* Members */ /* ------------------------------------------------------------------------ */ private: /// solid mechanics model SolidMechanicsModel * solid{nullptr}; /// phasefield model PhaseFieldModel * phase{nullptr}; Array * displacement{nullptr}; /// Array * displacement_increment{nullptr}; /// external forces array Array * external_force{nullptr}; }; } // namespace akantu #endif /* __AKANTU_COUPLER_SOLID_PHASEFIELD_HH__ */ diff --git a/src/model/phase_field/phase_field_model.cc b/src/model/phase_field/phase_field_model.cc index e35be4a8c..5a1988997 100644 --- a/src/model/phase_field/phase_field_model.cc +++ b/src/model/phase_field/phase_field_model.cc @@ -1,776 +1,776 @@ /** * @file phase_field_model.cc * * @author Mohit Pundir * * @date creation: Wed Aug 01 2018 * @date last modification: Wed Aug 01 2018 * * @brief Implementation of PhaseFieldModel class * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "phase_field_model.hh" #include "dumpable_inline_impl.hh" #include "element_synchronizer.hh" #include "fe_engine_template.hh" #include "generalized_trapezoidal.hh" #include "group_manager_inline_impl.hh" #include "integrator_gauss.hh" #include "mesh.hh" #include "parser.hh" #include "shape_lagrange.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_element_partition.hh" #include "dumper_elemental_field.hh" #include "dumper_internal_material_field.hh" #include "dumper_iohelper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ PhaseFieldModel::PhaseFieldModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id, const ModelType model_type) : Model(mesh, model_type, dim, id, memory_id), phasefield_index("phasefield index", id, memory_id), phasefield_local_numbering("phasefield local numbering", id, memory_id) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("PhaseFieldFEEngine", mesh, Model::spatial_dimension); #ifdef AKANTU_USE_IOHELPER this->mesh.registerDumper("phase_field", id, true); this->mesh.addDumpMesh(mesh, Model::spatial_dimension, _not_ghost, _ek_regular); #endif // AKANTU_USE_IOHELPER phasefield_selector = std::make_shared(phasefield_index); this->initDOFManager(); this->registerDataAccessor(*this); if (this->mesh.isDistributed()) { auto & synchronizer = this->mesh.getElementSynchronizer(); this->registerSynchronizer(synchronizer, SynchronizationTag::_pfm_damage); this->registerSynchronizer(synchronizer, SynchronizationTag::_pfm_driving); this->registerSynchronizer(synchronizer, SynchronizationTag::_pfm_history); this->registerSynchronizer(synchronizer, SynchronizationTag::_pfm_energy); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ PhaseFieldModel::~PhaseFieldModel() = default; /* -------------------------------------------------------------------------- */ MatrixType PhaseFieldModel::getMatrixType(const ID & matrix_id) { if (matrix_id == "K" or matrix_id == "M") { return _symmetric; } return _mt_not_defined; } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::initModel() { auto & fem = this->getFEEngine(); fem.initShapeFunctions(_not_ghost); fem.initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::initFullImpl(const ModelOptions & options) { phasefield_index.initialize(mesh, _element_kind = _ek_not_defined, _default_value = UInt(-1), _with_nb_element = true); phasefield_local_numbering.initialize(mesh, _element_kind = _ek_not_defined, _with_nb_element = true); Model::initFullImpl(options); // initialize the phasefields if (this->parser.getLastParsedFile() != "") { this->instantiatePhaseFields(); this->initPhaseFields(); } this->initBC(*this, *damage, *external_force); } /* -------------------------------------------------------------------------- */ PhaseField & PhaseFieldModel::registerNewPhaseField(const ParserSection & section) { std::string phase_name; std::string phase_type = section.getName(); std::string opt_param = section.getOption(); try { std::string tmp = section.getParameter("name"); phase_name = tmp; /** this can seam weird, but there is an ambiguous operator * overload that i couldn't solve. @todo remove the * weirdness of this code */ } catch (debug::Exception &) { AKANTU_ERROR("A phasefield of type \'" << phase_type << "\' in the input file has been defined without a name!"); } PhaseField & phase = this->registerNewPhaseField(phase_name, phase_type, opt_param); phase.parseSection(section); return phase; } /* -------------------------------------------------------------------------- */ PhaseField & PhaseFieldModel::registerNewPhaseField(const ID & phase_name, const ID & phase_type, const ID & opt_param) { AKANTU_DEBUG_ASSERT(phasefields_names_to_id.find(phase_name) == phasefields_names_to_id.end(), "A phasefield with this name '" << phase_name << "' has already been registered. " << "Please use unique names for phasefields"); UInt phase_count = phasefields.size(); phasefields_names_to_id[phase_name] = phase_count; std::stringstream sstr_phase; sstr_phase << this->id << ":" << phase_count << ":" << phase_type; ID mat_id = sstr_phase.str(); std::unique_ptr phase = PhaseFieldFactory::getInstance().allocate( phase_type, spatial_dimension, opt_param, *this, mat_id); phasefields.push_back(std::move(phase)); return *(phasefields.back()); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::instantiatePhaseFields() { ParserSection model_section; bool is_empty; std::tie(model_section, is_empty) = this->getParserSection(); if (not is_empty) { auto model_phasefields = model_section.getSubSections(ParserType::_phasefield); for (const auto & section : model_phasefields) { this->registerNewPhaseField(section); } } auto sub_sections = this->parser.getSubSections(ParserType::_phasefield); for (const auto & section : sub_sections) { this->registerNewPhaseField(section); } if (phasefields.empty()) AKANTU_EXCEPTION("No phasefields where instantiated for the model" << getID()); are_phasefields_instantiated = true; } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::initPhaseFields() { AKANTU_DEBUG_ASSERT(phasefields.size() != 0, "No phasefield to initialize !"); if (!are_phasefields_instantiated) instantiatePhaseFields(); this->assignPhaseFieldToElements(); for (auto & phasefield : phasefields) { /// init internals properties phasefield->initPhaseField(); } this->synchronize(SynchronizationTag::_smm_init_mat); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::assignPhaseFieldToElements( const ElementTypeMapArray * filter) { for_each_element( mesh, [&](auto && element) { UInt phase_index = (*phasefield_selector)(element); AKANTU_DEBUG_ASSERT( phase_index < phasefields.size(), "The phasefield selector returned an index that does not exists"); phasefield_index(element) = phase_index; }, _element_filter = filter, _ghost_type = _not_ghost); for_each_element(mesh, [&](auto && element) { auto phase_index = phasefield_index(element); auto index = phasefields[phase_index]->addElement(element); phasefield_local_numbering(element) = index; }, _element_filter = filter, _ghost_type = _not_ghost); // synchronize the element phasefield arrays this->synchronize(SynchronizationTag::_material_id); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } else { AKANTU_ERROR("Unknown Matrix ID for PhaseFieldModel : " << matrix_id); } } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::predictor() { // AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::corrector() { // AKANTU_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType) { DOFManager & dof_manager = this->getDOFManager(); this->allocNodalField(this->damage, 1, "damage"); this->allocNodalField(this->external_force, 1, "external_force"); this->allocNodalField(this->internal_force, 1, "internal_force"); this->allocNodalField(this->blocked_dofs, 1, "blocked_dofs"); this->allocNodalField(this->previous_damage, 1, "previous_damage"); this->allocNodalField(this->damage_increment, 1, "damage_increment"); if (!dof_manager.hasDOFs("damage")) { dof_manager.registerDOFs("damage", *this->damage, _dst_nodal); dof_manager.registerBlockedDOFs("damage", *this->blocked_dofs); dof_manager.registerDOFsIncrement("damage", *this->damage_increment); dof_manager.registerDOFsPrevious("damage", *this->previous_damage); } if (time_step_solver_type == TimeStepSolverType::_dynamic) { AKANTU_TO_IMPLEMENT(); } } /* -------------------------------------------------------------------------- */ FEEngine & PhaseFieldModel::getFEEngineBoundary(const ID & name) { return dynamic_cast(getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ std::tuple PhaseFieldModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _explicit_lumped_mass: { return std::make_tuple("explicit_lumped", TimeStepSolverType::_dynamic_lumped); } case _explicit_consistent_mass: { return std::make_tuple("explicit", TimeStepSolverType::_dynamic); } case _static: { return std::make_tuple("static", TimeStepSolverType::_static); } case _implicit_dynamic: { return std::make_tuple("implicit", TimeStepSolverType::_dynamic); } default: return std::make_tuple("unknown", TimeStepSolverType::_not_defined); } } /* -------------------------------------------------------------------------- */ ModelSolverOptions PhaseFieldModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_dynamic_lumped: { options.non_linear_solver_type = NonLinearSolverType::_lumped; options.integration_scheme_type["damage"] = IntegrationSchemeType::_central_difference; options.solution_type["damage"] = IntegrationScheme::_acceleration; break; } case TimeStepSolverType::_static: { options.non_linear_solver_type = NonLinearSolverType::_linear; options.integration_scheme_type["damage"] = IntegrationSchemeType::_pseudo_time; options.solution_type["damage"] = IntegrationScheme::_not_defined; break; } case TimeStepSolverType::_dynamic: { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["damage"] = IntegrationSchemeType::_backward_euler; options.solution_type["damage"] = IntegrationScheme::_damage; break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); } return options; } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::beforeSolveStep() { for (auto & phasefield : phasefields) { phasefield->beforeSolveStep(); } // compute the history of local elements //AKANTU_DEBUG_INFO("Compute phi history"); //this->computePhiHistoryOnQuadPoints(_not_ghost); // communicate the history //AKANTU_DEBUG_INFO("Send data for synchronization"); //this->asynchronousSynchronize(SynchronizationTag::_pfm_history); // finalize communications //AKANTU_DEBUG_INFO("Wait distant history"); //this->waitEndSynchronize(SynchronizationTag::_pfm_history); //this->computeDamageEnergyDensityOnQuadPoints(_not_ghost); // communicate the energy density //AKANTU_DEBUG_INFO("Send data for synchronization"); //this->asynchronousSynchronize(SynchronizationTag::_pfm_energy); // finalize communications //AKANTU_DEBUG_INFO("Wait distant energy density"); //this->waitEndSynchronize(SynchronizationTag::_pfm_energy); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::afterSolveStep(bool converged) { if (not converged) return ; for (auto && values : zip(*damage, *previous_damage)) { auto & dam = std::get<0>(values); auto & prev_dam = std::get<1>(values); dam -= prev_dam; //dam = std::min(1., 2 * dam - dam * dam); prev_dam = dam; } } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix"); if (!this->getDOFManager().hasMatrix("K")) { this->getDOFManager().getNewMatrix("K", getMatrixType("K")); } this->getDOFManager().zeroMatrix("K"); for (auto & phasefield : phasefields) { phasefield->assembleStiffnessMatrix(_not_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::assembleResidual() { AKANTU_DEBUG_IN(); this->assembleInternalForces(); this->getDOFManager().assembleToResidual("damage", *this->external_force, 1); this->getDOFManager().assembleToResidual("damage", *this->internal_force, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); this->internal_force->zero(); // compute the driving forces of local elements AKANTU_DEBUG_INFO("Compute local driving forces"); for (auto & phasefield : phasefields) { phasefield->computeAllDrivingForces(_not_ghost); } // communicate the driving forces AKANTU_DEBUG_INFO("Send data for residual assembly"); this->asynchronousSynchronize(SynchronizationTag::_pfm_driving); // assemble the forces due to local driving forces AKANTU_DEBUG_INFO("Assemble residual for local elements"); for (auto & phasefield : phasefields) { phasefield->assembleInternalForces(_not_ghost); } // finalize communications AKANTU_DEBUG_INFO("Wait distant driving forces"); this->waitEndSynchronize(SynchronizationTag::_pfm_driving); // assemble the residual due to ghost elements AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); //for (auto & phasefield : phasefields) { // phasefield->assembleInternalForces(_ghost); //} AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::assembleLumpedMatrix(const ID & /*matrix_id*/) {} /* -------------------------------------------------------------------------- */ void PhaseFieldModel::setTimeStep(Real time_step, const ID & solver_id) { Model::setTimeStep(time_step, solver_id); #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper("phase_field").setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ UInt PhaseFieldModel::getNbData(const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes_per_element = 0; for (const Element & el : elements) { nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); } switch (tag) { case SynchronizationTag::_pfm_damage: { size += nb_nodes_per_element * sizeof(Real); // damage break; } case SynchronizationTag::_pfm_driving: { size += getNbIntegrationPoints(elements) * sizeof(Real); break; } case SynchronizationTag::_pfm_history: { size += getNbIntegrationPoints(elements) * sizeof(Real); break; } case SynchronizationTag::_pfm_energy: { size += getNbIntegrationPoints(elements) * sizeof(Real); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ -void PhaseFieldModel::packData(CommunicationBuffer & buffer, - const Array & elements, - const SynchronizationTag & tag) const { +void PhaseFieldModel::packData(__attribute__((unused)) CommunicationBuffer & buffer, + __attribute__((unused)) const Array & elements, + __attribute__((unused)) const SynchronizationTag & tag) const { /*switch (tag) { case SynchronizationTag::_pfm_damage: { packNodalDataHelper(*damage, buffer, elements, mesh); break; } case SynchronizationTag::_pfm_driving: { packElementalDataHelper(driving_force_on_qpoints, buffer, elements, true, getFEEngine()); break; } case SynchronizationTag::_pfm_history: { packElementalDataHelper(phi_history_on_qpoints, buffer, elements, true, getFEEngine()); break; } case SynchronizationTag::_pfm_energy: { packElementalDataHelper(damage_energy_density_on_qpoints, buffer, elements, true, getFEEngine()); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } }*/ } /* -------------------------------------------------------------------------- */ -void PhaseFieldModel::unpackData(CommunicationBuffer & buffer, - const Array & elements, - const SynchronizationTag & tag) { +void PhaseFieldModel::unpackData(__attribute__((unused)) CommunicationBuffer & buffer, + __attribute__((unused)) const Array & elements, + __attribute__((unused)) const SynchronizationTag & tag) { /*switch (tag) { case SynchronizationTag::_pfm_damage: { unpackNodalDataHelper(*damage, buffer, elements, mesh); break; } case SynchronizationTag::_pfm_driving: { unpackElementalDataHelper(driving_force_on_qpoints, buffer, elements, true, getFEEngine()); break; } case SynchronizationTag::_pfm_history: { unpackElementalDataHelper(phi_history_on_qpoints, buffer, elements, true, getFEEngine()); break; } case SynchronizationTag::_pfm_energy: { unpackElementalDataHelper(damage_energy_density_on_qpoints, buffer, elements, true, getFEEngine()); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } }*/ } /* -------------------------------------------------------------------------- */ UInt PhaseFieldModel::getNbData(const Array & indexes, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes = indexes.size(); switch (tag) { case SynchronizationTag::_pfm_damage: { size += nb_nodes * sizeof(Real); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::packData(CommunicationBuffer & buffer, const Array & indexes, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); for (auto index : indexes) { switch (tag) { case SynchronizationTag::_pfm_damage: { buffer << (*damage)(index); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::unpackData(CommunicationBuffer & buffer, const Array & indexes, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); for (auto index : indexes) { switch (tag) { case SynchronizationTag::_pfm_damage: { buffer >> (*damage)(index); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER std::shared_ptr PhaseFieldModel::createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool) { std::map *> uint_nodal_fields; uint_nodal_fields["blocked_dofs"] = blocked_dofs; return mesh.createNodalField(uint_nodal_fields[field_name], group_name); std::shared_ptr field; return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr PhaseFieldModel::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool) { 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["damage"] = damage; real_nodal_fields["external_force"] = external_force; real_nodal_fields["internal_force"] = internal_force; return mesh.createNodalField(real_nodal_fields[field_name], group_name); std::shared_ptr field; return field; } /* -------------------------------------------------------------------------- */ std::shared_ptr PhaseFieldModel::createElementalField(const std::string & field_name, const std::string & group_name, bool, UInt, ElementKind element_kind) { if (field_name == "partitions") { return mesh.createElementalField( mesh.getConnectivities(), group_name, this->spatial_dimension, element_kind); } /*else if (field_name == "damage_gradient") { ElementTypeMap nb_data_per_elem = this->mesh.getNbDataPerElem(damage_gradient); return mesh.createElementalField( damage_gradient, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); } else if (field_name == "damage_energy") { ElementTypeMap nb_data_per_elem = this->mesh.getNbDataPerElem(damage_energy_on_qpoints); return mesh.createElementalField( damage_energy_on_qpoints, group_name, this->spatial_dimension, element_kind, nb_data_per_elem); }*/ std::shared_ptr field; return field; } /* -------------------------------------------------------------------------- */ #else /* -------------------------------------------------------------------------- */ std::shared_ptr PhaseFieldModel::createElementalField(const std::string &, const std::string &, bool, const UInt &, ElementKind) { return nullptr; } /* -------------------------------------------------------------------------- */ std::shared_ptr PhaseFieldModel::createNodalFieldReal(const std::string &, const std::string &, bool) { return nullptr; } /* -------------------------------------------------------------------------- */ std::shared_ptr PhaseFieldModel::createNodalFieldBool(const std::string &, const std::string &, bool) { return nullptr; } #endif /* -------------------------------------------------------------------------- */ void PhaseFieldModel::dump(const std::string & dumper_name) { mesh.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::dump(const std::string & dumper_name, UInt step) { mesh.dump(dumper_name, step); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::dump(const std::string & dumper_name, Real time, UInt step) { mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::dump() { mesh.dump(); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::dump(UInt step) { mesh.dump(step); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::dump(Real time, UInt step) { mesh.dump(time, step); } /* -------------------------------------------------------------------------- */ void PhaseFieldModel::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "Phase Field Model [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + spatial dimension : " << Model::spatial_dimension << std::endl; stream << space << " + fem [" << std::endl; getFEEngine().printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + nodals information [" << std::endl; damage->printself(stream, indent + 2); external_force->printself(stream, indent + 2); internal_force->printself(stream, indent + 2); blocked_dofs->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + phasefield information [" << std::endl; stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } } // namespace akantu diff --git a/src/model/phase_field/phasefield.cc b/src/model/phase_field/phasefield.cc index 529c18679..538541cf2 100644 --- a/src/model/phase_field/phasefield.cc +++ b/src/model/phase_field/phasefield.cc @@ -1,302 +1,302 @@ /** * @file pahsefield.cc * * @author Mohit Pundir * * @date creation: Mon Mar 2 2020 * @date last modification: Mon Mar 2 2020 * * @brief Implementation of the common part of the phasefield class * * @section LICENSE * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "phasefield.hh" #include "phase_field_model.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ PhaseField::PhaseField(PhaseFieldModel & model, const ID & id) : Memory(id, model.getMemoryID()), Parsable(ParserType::_phasefield, id), fem(model.getFEEngine()), name(""), model(model), spatial_dimension(this->model.getSpatialDimension()), element_filter("element_filter", id, this->memory_id), damage("damage", *this), phi("phi", *this), strain("strain", *this), driving_force("driving_force", *this), damage_energy("damage_energy", *this), damage_energy_density("damage_energy_density", *this) { AKANTU_DEBUG_IN(); /// for each connectivity types allocate the element filer array of the /// material element_filter.initialize(model.getMesh(), _spatial_dimension = spatial_dimension, _element_kind = _ek_regular); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ PhaseField::PhaseField(PhaseFieldModel & model, UInt dim, const Mesh & mesh, FEEngine & fe_engine, const ID & id) : Memory(id, model.getMemoryID()), Parsable(ParserType::_phasefield, id), fem(fe_engine), name(""), model(model), spatial_dimension(this->model.getSpatialDimension()), element_filter("element_filter", id, this->memory_id), damage("damage", *this, dim, fe_engine, this->element_filter), phi("phi", *this, dim, fe_engine, this->element_filter), strain("strain", *this, dim, fe_engine, this->element_filter), driving_force("driving_force", *this, dim, fe_engine, this->element_filter), damage_energy("damage_energy", *this, dim, fe_engine, this->element_filter), damage_energy_density("damage_energy_density", *this, dim, fe_engine, this->element_filter){ AKANTU_DEBUG_IN(); /// for each connectivity types allocate the element filer array of the /// material element_filter.initialize(mesh, _spatial_dimension = spatial_dimension, _element_kind = _ek_regular); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ PhaseField::~PhaseField() = default; /* -------------------------------------------------------------------------- */ void PhaseField::initialize() { registerParam("name", name, std::string(), _pat_parsable | _pat_readable); registerParam("l0", l0, Real(0.), _pat_parsable | _pat_readable, "length scale parameter"); registerParam("gc", g_c, _pat_parsable | _pat_readable, "critical local fracture energy density"); registerParam("E", E, _pat_parsable | _pat_readable, "Young's modulus"); registerParam("nu", nu, _pat_parsable | _pat_readable, "Poisson ratio"); damage.initialize(1); phi.initialize(1); driving_force.initialize(1); strain.initialize(spatial_dimension * spatial_dimension); damage_energy_density.initialize(1); damage_energy.initialize(spatial_dimension * spatial_dimension); } /* -------------------------------------------------------------------------- */ void PhaseField::initPhaseField() { AKANTU_DEBUG_IN(); this->phi.initializeHistory(); this->resizeInternals(); updateInternalParameters(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseField::resizeInternals() { AKANTU_DEBUG_IN(); for (auto it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) it->second->resize(); for (auto it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) it->second->resize(); for (auto it = internal_vectors_bool.begin(); it != internal_vectors_bool.end(); ++it) it->second->resize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseField::updateInternalParameters() { this->lambda = this->nu * this->E / ((1 + this->nu) * (1 - 2 * this->nu)); this->mu = this->E / (2 * (1 + this->nu)); } /* -------------------------------------------------------------------------- */ void PhaseField::computeAllDrivingForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); for (const auto & type : element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_filter = element_filter(type, ghost_type); if (elem_filter.size() == 0) continue; computeDrivingForce(type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseField::assembleInternalForces(GhostType ghost_type) { AKANTU_DEBUG_IN(); Array & internal_force = model.getInternalForce(); for (auto type : element_filter.elementTypes(_ghost_type = ghost_type)) { Array & elem_filter = element_filter(type, ghost_type); if (elem_filter.size() == 0) continue; UInt nb_element = elem_filter.size(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); Array nt_driving_force(nb_quadrature_points, nb_nodes_per_element); fem.computeNtb(driving_force(type, ghost_type), nt_driving_force, type, ghost_type, elem_filter); Array int_nt_driving_force(nb_element, nb_nodes_per_element); fem.integrate(nt_driving_force, int_nt_driving_force, nb_nodes_per_element, type, ghost_type, elem_filter); model.getDOFManager().assembleElementalArrayLocalArray( int_nt_driving_force, internal_force, type, ghost_type, 1, elem_filter); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseField::assembleStiffnessMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix"); for (auto type : element_filter.elementTypes(spatial_dimension, ghost_type)) { Array & elem_filter = element_filter(type, ghost_type); if (elem_filter.size() == 0) { AKANTU_DEBUG_OUT(); return; } auto nb_element = elem_filter.size(); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_quadrature_points = fem.getNbIntegrationPoints(type, ghost_type); auto nt_b_n = std::make_unique>( nb_element * nb_quadrature_points, nb_nodes_per_element * nb_nodes_per_element, "N^t*b*N"); auto bt_d_b = std::make_unique>( nb_element * nb_quadrature_points, nb_nodes_per_element * nb_nodes_per_element, "B^t*D*B"); - // damage_energy_on_qpoints = gc/l0 + phi = scalar + // damage_energy_density_on_qpoints = gc/l0 + phi = scalar auto & damage_energy_density_vect = damage_energy_density(type, ghost_type); // damage_energy_on_qpoints = gc*l0 = scalar auto & damage_energy_vect = damage_energy(type, ghost_type); fem.computeBtDB(damage_energy_vect, *bt_d_b, 2, type, ghost_type, elem_filter); fem.computeNtbN(damage_energy_density_vect, *nt_b_n, 2, type, ghost_type, elem_filter); /// compute @f$ K_{\grad d} = \int_e \mathbf{N}^t * \mathbf{w} * /// \mathbf{N}@f$ auto K_n = std::make_unique>( nb_element, nb_nodes_per_element * nb_nodes_per_element, "K_n"); fem.integrate(*nt_b_n, *K_n, nb_nodes_per_element * nb_nodes_per_element, type, ghost_type, elem_filter); model.getDOFManager().assembleElementalMatricesToMatrix( "K", "damage", *K_n, type, _not_ghost, _symmetric, elem_filter); /// compute @f$ K_{\grad d} = \int_e \mathbf{B}^t * \mathbf{W} * /// \mathbf{B}@f$ auto K_b = std::make_unique>( nb_element, nb_nodes_per_element * nb_nodes_per_element, "K_b"); fem.integrate(*bt_d_b, *K_b, nb_nodes_per_element * nb_nodes_per_element, type, ghost_type, elem_filter); model.getDOFManager().assembleElementalMatricesToMatrix( "K", "damage", *K_b, type, _not_ghost, _symmetric, elem_filter); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseField::beforeSolveStep() { this->savePreviousState(); } /* -------------------------------------------------------------------------- */ void PhaseField::afterSolveStep() { } /* -------------------------------------------------------------------------- */ void PhaseField::savePreviousState() { AKANTU_DEBUG_IN(); for (auto pair : internal_vectors_real) if (pair.second->hasHistory()) pair.second->saveCurrentValues(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void PhaseField::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); std::string type = getID().substr(getID().find_last_of(':') + 1); stream << space << "PhaseField Material " << type << " [" << std::endl; Parsable::printself(stream, indent); stream << space << "]" << std::endl; } } diff --git a/src/model/phase_field/phasefield_inline_impl.cc b/src/model/phase_field/phasefield_inline_impl.cc index af7d2d54a..960c64e22 100644 --- a/src/model/phase_field/phasefield_inline_impl.cc +++ b/src/model/phase_field/phasefield_inline_impl.cc @@ -1,147 +1,147 @@ /* -------------------------------------------------------------------------- */ #include "phase_field_model.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_PHASEFIELD_INLINE_IMPL_CC__ #define __AKANTU_PHASEFIELD_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ inline UInt PhaseField::addElement(const ElementType & type, UInt element, const GhostType & ghost_type) { Array & el_filter = this->element_filter(type, ghost_type); el_filter.push_back(element); return el_filter.size() - 1; } /* -------------------------------------------------------------------------- */ inline UInt PhaseField::addElement(const Element & element) { return this->addElement(element.type, element.element, element.ghost_type); } /* -------------------------------------------------------------------------- */ template <> inline void PhaseField::registerInternal(InternalPhaseField & vect) { internal_vectors_real[vect.getID()] = &vect; } template <> inline void PhaseField::registerInternal(InternalPhaseField & vect) { internal_vectors_uint[vect.getID()] = &vect; } template <> inline void PhaseField::registerInternal(InternalPhaseField & vect) { internal_vectors_bool[vect.getID()] = &vect; } /* -------------------------------------------------------------------------- */ template <> inline void PhaseField::unregisterInternal(InternalPhaseField & vect) { internal_vectors_real.erase(vect.getID()); } template <> inline void PhaseField::unregisterInternal(InternalPhaseField & vect) { internal_vectors_uint.erase(vect.getID()); } template <> inline void PhaseField::unregisterInternal(InternalPhaseField & vect) { internal_vectors_bool.erase(vect.getID()); } /* -------------------------------------------------------------------------- */ template inline bool PhaseField::isInternal(__attribute__((unused)) const ID & id, __attribute__((unused)) const ElementKind & element_kind) const { AKANTU_TO_IMPLEMENT(); } template <> inline bool PhaseField::isInternal(const ID & id, const ElementKind & element_kind) const { auto internal_array = internal_vectors_real.find(this->getID() + ":" + id); if (internal_array == internal_vectors_real.end() || internal_array->second->getElementKind() != element_kind) return false; return true; } /* -------------------------------------------------------------------------- */ -inline UInt PhaseField::getNbData(const Array & elements, - const SynchronizationTag & tag) const { +inline UInt PhaseField::getNbData(__attribute__((unused)) const Array & elements, + __attribute__((unused)) const SynchronizationTag & tag) const { /*if (tag == SynchronizationTag::_smm_stress) { return (this->isFiniteDeformation() ? 3 : 1) * spatial_dimension * spatial_dimension * sizeof(Real) * this->getModel().getNbIntegrationPoints(elements); }*/ return 0; } /* -------------------------------------------------------------------------- */ -inline void PhaseField::packData(CommunicationBuffer & buffer, - const Array & elements, - const SynchronizationTag & tag) const { +inline void PhaseField::packData(__attribute__((unused)) CommunicationBuffer & buffer, + __attribute__((unused)) const Array & elements, + __attribute__((unused)) const SynchronizationTag & tag) const { /*if (tag == SynchronizationTag::_smm_stress) { if (this->isFiniteDeformation()) { packElementDataHelper(piola_kirchhoff_2, buffer, elements); packElementDataHelper(gradu, buffer, elements); } packElementDataHelper(stress, buffer, elements); }*/ } /* -------------------------------------------------------------------------- */ -inline void PhaseField::unpackData(CommunicationBuffer & buffer, - const Array & elements, - const SynchronizationTag & tag) { +inline void PhaseField::unpackData(__attribute__((unused)) CommunicationBuffer & buffer, + __attribute__((unused)) const Array & elements, + __attribute__((unused)) const SynchronizationTag & tag) { /*if (tag == SynchronizationTag::_smm_stress) { if (this->isFiniteDeformation()) { unpackElementDataHelper(piola_kirchhoff_2, buffer, elements); unpackElementDataHelper(gradu, buffer, elements); } unpackElementDataHelper(stress, buffer, elements); }*/ } /* -------------------------------------------------------------------------- */ inline const Parameter & PhaseField::getParam(const ID & param) const { try { return get(param); } catch (...) { AKANTU_EXCEPTION("No parameter " << param << " in the material " << getID()); } } /* -------------------------------------------------------------------------- */ template inline void PhaseField::packElementDataHelper( const ElementTypeMapArray & data_to_pack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id) const { DataAccessor::packElementalDataHelper(data_to_pack, buffer, elements, true, model.getFEEngine(fem_id)); } /* -------------------------------------------------------------------------- */ template inline void PhaseField::unpackElementDataHelper( ElementTypeMapArray & data_to_unpack, CommunicationBuffer & buffer, const Array & elements, const ID & fem_id) { DataAccessor::unpackElementalDataHelper(data_to_unpack, buffer, elements, true, model.getFEEngine(fem_id)); } } #endif diff --git a/test/test_model/test_phase_field_model/test_phase_solid_coupling.cc b/test/test_model/test_phase_field_model/test_phase_solid_coupling.cc index 03cc09ad6..ad53071ad 100644 --- a/test/test_model/test_phase_field_model/test_phase_solid_coupling.cc +++ b/test/test_model/test_phase_field_model/test_phase_solid_coupling.cc @@ -1,280 +1,272 @@ /** * @file test_phase_field_coupling.cc * * @author Mohit Pundir * * @date creation: Thu Feb 25 2021 * * @brief test of the class PhaseFieldModel on the 2d square * * @section LICENSE * * Copyright (©) 2015 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_common.hh" #include "non_linear_solver.hh" #include "solid_mechanics_model.hh" #include "phase_field_model.hh" #include "material.hh" #include "material_phasefield.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; const UInt spatial_dimension = 2; /* -------------------------------------------------------------------------- */ void applyDisplacement(SolidMechanicsModel &, Real &); void computeStrainOnQuadPoints(SolidMechanicsModel &, PhaseFieldModel &, const GhostType &); void computeDamageOnQuadPoints(SolidMechanicsModel &, PhaseFieldModel &, const GhostType &); void gradUToEpsilon(const Matrix &, Matrix &); /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { std::ofstream os("data.csv"); os << "#strain stress damage analytical_sigma analytical_damage" << std::endl; initialize("material_coupling.dat", argc, argv); Mesh mesh(spatial_dimension); mesh.read("test_one_element.msh"); SolidMechanicsModel model(mesh); model.initFull(_analysis_method = _static); PhaseFieldModel phase(mesh); auto && selector = std::make_shared>( "physical_names", phase); phase.setPhaseFieldSelector(selector); phase.initFull(_analysis_method = _static); + model.setBaseName("phase_solid"); model.addDumpField("stress"); model.addDumpField("grad_u"); model.addDumpFieldVector("displacement"); model.addDumpField("damage"); model.dump(); UInt nbSteps = 1000; Real increment = 1e-4; auto & stress = model.getMaterial(0).getArray("stress", _quadrangle_4); auto & damage = model.getMaterial(0).getArray("damage", _quadrangle_4); Real analytical_damage{0.}; Real analytical_sigma{0.}; - auto & mat = model.getMaterial(0); auto & phasefield = phase.getPhaseField(0); const Real E = phasefield.getParam("E"); const Real nu = phasefield.getParam("nu"); Real c22 = E*(1-nu)/((1+nu)*(1-2*nu)); const Real gc = phasefield.getParam("gc"); const Real l0 = phasefield.getParam("l0"); - Real error_stress{0.}; Real error_damage{0.}; for (UInt s = 0; s < nbSteps; ++s) { Real axial_strain = increment * s; applyDisplacement(model, axial_strain); model.solveStep(); computeStrainOnQuadPoints(model, phase, _not_ghost); phase.solveStep(); computeDamageOnQuadPoints(model, phase, _not_ghost); model.assembleInternalForces(); analytical_damage = axial_strain*axial_strain*c22/(gc/l0 + axial_strain*axial_strain*c22); analytical_sigma = c22*axial_strain*(1-analytical_damage)*(1-analytical_damage); - - error_stress = std::abs(analytical_sigma - stress(0, 3))/analytical_sigma; error_damage = std::abs(analytical_damage - damage(0))/analytical_damage; if (error_damage > 0.01) { return EXIT_FAILURE; } os << axial_strain << " " << stress(0, 3) << " " << damage(0) << " " << analytical_sigma << " " << analytical_damage << std::endl; model.dump(); } os.close(); finalize(); return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */ void applyDisplacement(SolidMechanicsModel & model, Real & increment) { auto & displacement = model.getDisplacement(); auto & positions = model.getMesh().getNodes(); auto & blocked_dofs = model.getBlockedDOFs(); for (UInt n = 0; n < model.getMesh().getNbNodes(); ++n) { if (positions(n, 1) == -0.5) { displacement(n, 0) = 0; displacement(n, 1) = 0; blocked_dofs(n, 0) = true; blocked_dofs(n ,1) = true; } else { displacement(n, 0) = 0; displacement(n, 1) = increment; blocked_dofs(n, 0) = true; blocked_dofs(n ,1) = true; } } } /* -------------------------------------------------------------------------- */ void computeStrainOnQuadPoints(SolidMechanicsModel & solid, PhaseFieldModel & phase, const GhostType & ghost_type) { auto & mesh = solid.getMesh(); auto nb_materials = solid.getNbMaterials(); auto nb_phasefields = phase.getNbPhaseFields(); AKANTU_DEBUG_ASSERT(nb_phasefields == nb_materials, "The number of phasefields and materials should be equal" ); for(auto index : arange(nb_materials)) { auto & material = solid.getMaterial(index); for(auto index2 : arange(nb_phasefields)) { auto & phasefield = phase.getPhaseField(index2); if(phasefield.getName() == material.getName()){ auto & strain_on_qpoints = phasefield.getStrain(); auto & gradu_on_qpoints = material.getGradU(); for (auto & type: mesh.elementTypes(spatial_dimension, ghost_type)) { auto & strain_on_qpoints_vect = strain_on_qpoints(type, ghost_type); auto & gradu_on_qpoints_vect = gradu_on_qpoints(type, ghost_type); for (auto && values: zip(make_view(strain_on_qpoints_vect, spatial_dimension, spatial_dimension), make_view(gradu_on_qpoints_vect, spatial_dimension, spatial_dimension))) { auto & strain = std::get<0>(values); auto & grad_u = std::get<1>(values); gradUToEpsilon(grad_u, strain); } } break; } } } } /* -------------------------------------------------------------------------- */ void computeDamageOnQuadPoints(SolidMechanicsModel & solid, PhaseFieldModel & phase, const GhostType & ghost_type) { auto & fem = phase.getFEEngine(); auto & mesh = phase.getMesh(); auto nb_materials = solid.getNbMaterials(); auto nb_phasefields = phase.getNbPhaseFields(); AKANTU_DEBUG_ASSERT(nb_phasefields == nb_materials, "The number of phasefields and materials should be equal" ); for(auto index : arange(nb_materials)) { auto & material = solid.getMaterial(index); for(auto index2 : arange(nb_phasefields)) { auto & phasefield = phase.getPhaseField(index2); if(phasefield.getName() == material.getName()){ switch (spatial_dimension) { case 1: { auto & mat = static_cast &>(material); auto & solid_damage = mat.getDamage(); - auto & phase_damage = phasefield.getDamage(); for (auto & type: mesh.elementTypes(spatial_dimension, ghost_type)) { auto & damage_on_qpoints_vect = solid_damage(type, ghost_type); - auto & phase_damage_on_qpoints_vect = phase_damage(type, ghost_type); fem.interpolateOnIntegrationPoints(phase.getDamage(), damage_on_qpoints_vect, 1, type, ghost_type); } break; } case 2: { auto & mat = static_cast &>(material); auto & solid_damage = mat.getDamage(); - auto & phase_damage = phasefield.getDamage(); for (auto & type: mesh.elementTypes(spatial_dimension, ghost_type)) { auto & damage_on_qpoints_vect = solid_damage(type, ghost_type); - auto & phase_damage_on_qpoints_vect = phase_damage(type, ghost_type); fem.interpolateOnIntegrationPoints(phase.getDamage(), damage_on_qpoints_vect, 1, type, ghost_type); } break; } default: - auto & mat = static_cast &>(material); break; } } } } } /* -------------------------------------------------------------------------- */ void gradUToEpsilon(const Matrix & grad_u, Matrix & epsilon) { for (UInt i=0; i < spatial_dimension; ++i) { for (UInt j = 0; j < spatial_dimension; ++j) epsilon(i, j) = 0.5 * (grad_u(i, j) + grad_u(j, i)); } } diff --git a/test/test_model/test_phase_field_model/test_phase_solid_explicit.cc b/test/test_model/test_phase_field_model/test_phase_solid_explicit.cc index f9fc9eef9..b053f772c 100644 --- a/test/test_model/test_phase_field_model/test_phase_solid_explicit.cc +++ b/test/test_model/test_phase_field_model/test_phase_solid_explicit.cc @@ -1,161 +1,160 @@ /** * @file test_phase_field_explicit.cc * * @author Mohit Pundir * * @date creation: Thu Feb 28 2021 * * @brief test of the class PhaseFieldModel on the 2d square * * @section LICENSE * * Copyright (©) 2015 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_common.hh" #include "non_linear_solver.hh" #include "coupler_solid_phasefield.hh" #include "solid_mechanics_model.hh" #include "phase_field_model.hh" #include "material.hh" #include "material_phasefield.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; const UInt spatial_dimension = 2; /* -------------------------------------------------------------------------- */ void applyDisplacement(SolidMechanicsModel &, Real &); /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { std::ofstream os("data-explicit.csv"); os << "#strain stress damage analytical_sigma analytical_damage error_stress error_damage" << std::endl; initialize("material_coupling.dat", argc, argv); Mesh mesh(spatial_dimension); mesh.read("test_one_element.msh"); CouplerSolidPhaseField coupler(mesh); auto & model = coupler.getSolidMechanicsModel(); auto & phase = coupler.getPhaseFieldModel(); model.initFull(_analysis_method = _explicit_lumped_mass); Real time_factor = 0.8; Real stable_time_step = model.getStableTimeStep() * time_factor; model.setTimeStep(stable_time_step); auto && selector = std::make_shared>( "physical_names", phase); phase.setPhaseFieldSelector(selector); phase.initFull(_analysis_method = _static); model.setBaseName("phase_solid"); model.addDumpField("stress"); model.addDumpField("grad_u"); model.addDumpFieldVector("displacement"); model.addDumpField("damage"); model.dump(); UInt nbSteps = 1000; Real increment = 1e-4; auto & stress = model.getMaterial(0).getArray("stress", _quadrangle_4); auto & damage = model.getMaterial(0).getArray("damage", _quadrangle_4); Real analytical_damage{0.}; Real analytical_sigma{0.}; - auto & mat = model.getMaterial(0); auto & phasefield = phase.getPhaseField(0); const Real E = phasefield.getParam("E"); const Real nu = phasefield.getParam("nu"); Real c22 = E*(1-nu)/((1+nu)*(1-2*nu)); const Real gc = phasefield.getParam("gc"); const Real l0 = phasefield.getParam("l0"); Real error_stress{0.}; Real error_damage{0.}; for (UInt s = 0; s < nbSteps; ++s) { Real axial_strain = increment * s; applyDisplacement(model, axial_strain); - coupler.solve(); + coupler.solve("explicit_lumped", "static"); analytical_damage = axial_strain*axial_strain*c22/(gc/l0 + axial_strain*axial_strain*c22); analytical_sigma = c22*axial_strain*(1-analytical_damage)*(1-analytical_damage); error_stress = std::abs(analytical_sigma - stress(0, 3))/analytical_sigma; error_damage = std::abs(analytical_damage - damage(0))/analytical_damage; if (error_damage > 0.01) { return EXIT_FAILURE; } os << axial_strain << " " << stress(0, 3) << " " << damage(0) << " " << analytical_sigma << " " << analytical_damage << " " << error_stress << " " << error_damage << std::endl; model.dump(); } os.close(); finalize(); return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */ void applyDisplacement(SolidMechanicsModel & model, Real & increment) { auto & displacement = model.getDisplacement(); auto & positions = model.getMesh().getNodes(); auto & blocked_dofs = model.getBlockedDOFs(); for (UInt n = 0; n < model.getMesh().getNbNodes(); ++n) { if (positions(n, 1) == -0.5) { displacement(n, 0) = 0; displacement(n, 1) = 0; blocked_dofs(n, 0) = true; blocked_dofs(n ,1) = true; } else { displacement(n, 0) = 0; displacement(n, 1) = increment; blocked_dofs(n, 0) = true; blocked_dofs(n ,1) = true; } } }