diff --git a/python/py_model.cc b/python/py_model.cc index 32682d8a3..7f1c35e28 100644 --- a/python/py_model.cc +++ b/python/py_model.cc @@ -1,80 +1,81 @@ /* -------------------------------------------------------------------------- */ #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("getDOFManager", &Model::getDOFManager, py::return_value_policy::reference); } } // namespace akantu diff --git a/python/py_parser.cc b/python/py_parser.cc index bb0720ec0..2c8fc9dbf 100644 --- a/python/py_parser.cc +++ b/python/py_parser.cc @@ -1,69 +1,70 @@ /* -------------------------------------------------------------------------- */ #include "py_aka_array.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace py = pybind11; /* -------------------------------------------------------------------------- */ namespace akantu { std::map> map_params; void register_parser(py::module & mod) { py::enum_(mod, "ParameterAccessType", py::arithmetic()) .value("_pat_internal", _pat_internal) .value("_pat_writable", _pat_writable) .value("_pat_readable", _pat_readable) .value("_pat_modifiable", _pat_modifiable) .value("_pat_parsable", _pat_parsable) .value("_pat_parsmod", _pat_parsmod) .export_values(); py::class_(mod, "ParameterRegistry", py::multiple_inheritance()) .def("registerParamReal", [](ParameterRegistry & self, const std::string & name, UInt type, const std::string & description) { Real * p = new Real; map_params[&self][name] = p; self.registerParam(name, *p, ParameterAccessType(type), description); }) .def("registerParamReal", [](ParameterRegistry & self, const Real & _default, const std::string & name, UInt type, const std::string & description) { Real * p = new Real; map_params[&self][name] = p; self.registerParam(name, *p, _default, ParameterAccessType(type), description); }) .def("getReal", [](ParameterRegistry & self, const std::string & name) { return Real(self.get(name)); }) .def("getMatrix", [](ParameterRegistry & self, const std::string & name) { const Matrix & res = static_cast &>(self.get(name)); return res; }, py::return_value_policy::copy); py::class_(mod, "Parsable", py::multiple_inheritance()) .def(py::init()); mod.def("parseInput", [](const std::string & input_file) { getStaticParser().parse(input_file); }, "Parse an Akantu input file"); + } } // namespace akantu diff --git a/python/py_phase_field_model.cc b/python/py_phase_field_model.cc index 91cd30c42..b06c9609a 100644 --- a/python/py_phase_field_model.cc +++ b/python/py_phase_field_model.cc @@ -1,121 +1,122 @@ /* -------------------------------------------------------------------------- */ #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) { self.solve(); }) + .def("solve", [](CouplerSolidPhaseField & self, const ID & solver_id) { + self.solve(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 2972ed120..8c6761b74 100644 --- a/src/model/model_couplers/coupler_solid_phasefield.cc +++ b/src/model/model_couplers/coupler_solid_phasefield.cc @@ -1,704 +1,704 @@ /** * @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); 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() { +void CouplerSolidPhaseField::solve(const ID & solver_id) { - solid->solveStep(); + solid->solveStep(solver_id); this->computeStrainOnQuadPoints(_not_ghost); phase->solveStep(); 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, 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, 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 a6468d8c0..01f3aae09 100644 --- a/src/model/model_couplers/coupler_solid_phasefield.hh +++ b/src/model/model_couplers/coupler_solid_phasefield.hh @@ -1,289 +1,289 @@ /** * @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(); + void solve(const ID & 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 { 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 {} /* ------------------------------------------------------------------------ */ /* 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__ */