diff --git a/src/model/common/dof_manager/dof_manager_petsc.cc b/src/model/common/dof_manager/dof_manager_petsc.cc index 6fd1c286d..6cb7fb798 100644 --- a/src/model/common/dof_manager/dof_manager_petsc.cc +++ b/src/model/common/dof_manager/dof_manager_petsc.cc @@ -1,316 +1,316 @@ /** * @file dof_manager_petsc.cc * * @author Nicolas Richart * * @date creation: Wed Oct 07 2015 * @date last modification: Tue Feb 20 2018 * * @brief DOFManaterPETSc is the PETSc implementation of the DOFManager * * @section LICENSE * * Copyright (©) 2015-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 "dof_manager_petsc.hh" #include "aka_iterators.hh" #include "communicator.hh" #include "cppargparse.hh" #include "non_linear_solver_petsc.hh" #include "solver_vector_petsc.hh" #include "sparse_matrix_petsc.hh" #include "time_step_solver_default.hh" #if defined(AKANTU_USE_MPI) #include "mpi_communicator_data.hh" #endif /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace akantu { class PETScSingleton { private: PETScSingleton() { PETSc_call(PetscInitialized, &is_initialized); if (not is_initialized) { cppargparse::ArgumentParser & argparser = getStaticArgumentParser(); int & argc = argparser.getArgC(); char **& argv = argparser.getArgV(); PETSc_call(PetscInitialize, &argc, &argv, NULL, NULL); PETSc_call( PetscPopErrorHandler); // remove the default PETSc signal handler PETSc_call(PetscPushErrorHandler, PetscIgnoreErrorHandler, NULL); } } public: PETScSingleton(const PETScSingleton &) = delete; PETScSingleton & operator=(const PETScSingleton &) = delete; ~PETScSingleton() { if (not is_initialized) { PetscFinalize(); } } static PETScSingleton & getInstance() { static PETScSingleton instance; return instance; } private: PetscBool is_initialized; }; /* -------------------------------------------------------------------------- */ DOFManagerPETSc::DOFDataPETSc::DOFDataPETSc(const ID & dof_id) : DOFData(dof_id) {} /* -------------------------------------------------------------------------- */ DOFManagerPETSc::DOFManagerPETSc(const ID & id, const MemoryID & memory_id) : DOFManager(id, memory_id) { init(); } /* -------------------------------------------------------------------------- */ DOFManagerPETSc::DOFManagerPETSc(Mesh & mesh, const ID & id, const MemoryID & memory_id) : DOFManager(mesh, id, memory_id) { init(); } /* -------------------------------------------------------------------------- */ void DOFManagerPETSc::init() { // check if the akantu types and PETSc one are consistant static_assert(sizeof(Int) == sizeof(PetscInt), "The integer type of Akantu does not match the one from PETSc"); static_assert(sizeof(Real) == sizeof(PetscReal), "The integer type of Akantu does not match the one from PETSc"); #if defined(AKANTU_USE_MPI) const auto & mpi_data = aka::as_type(communicator.getCommunicatorData()); MPI_Comm mpi_comm = mpi_data.getMPICommunicator(); this->mpi_communicator = mpi_comm; #else this->mpi_communicator = PETSC_COMM_SELF; #endif PETScSingleton & instance [[gnu::unused]] = PETScSingleton::getInstance(); } /* -------------------------------------------------------------------------- */ DOFManagerPETSc::~DOFManagerPETSc() { // if (is_ltog_map) // PETSc_call(ISLocalToGlobalMappingDestroy, &is_ltog_map); } /* -------------------------------------------------------------------------- */ auto DOFManagerPETSc::getNewDOFData(const ID & dof_id) -> std::unique_ptr { return std::make_unique(dof_id); } /* -------------------------------------------------------------------------- */ std::tuple DOFManagerPETSc::registerDOFsInternal(const ID & dof_id, Array & dofs_array) { dofs_ids.push_back(dof_id); auto ret = DOFManager::registerDOFsInternal(dof_id, dofs_array); UInt nb_dofs, nb_pure_local_dofs; std::tie(nb_dofs, nb_pure_local_dofs, std::ignore) = ret; auto && vector = std::make_unique(*this, id + ":solution"); auto x = vector->getVec(); PETSc_call(VecGetLocalToGlobalMapping, x, &is_ltog_map); // redoing the indexes based on the petsc numbering for (auto & dof_id : dofs_ids) { auto & dof_data = this->getDOFDataTyped(dof_id); Array gidx(dof_data.local_equation_number.size()); for (auto && data : zip(dof_data.local_equation_number, gidx)) { std::get<1>(data) = localToGlobalEquationNumber(std::get<0>(data)); } auto & lidx = dof_data.local_equation_number_petsc; if (is_ltog_map) { lidx.resize(gidx.size()); PetscInt n; PETSc_call(ISGlobalToLocalMappingApply, is_ltog_map, IS_GTOLM_MASK, gidx.size(), gidx.storage(), &n, lidx.storage()); } } residual = std::make_unique(*vector, id + ":residual"); - data_cache = std::make_unique(*vector, id + ":residual"); + data_cache = std::make_unique(*vector, id + ":data_cache"); solution = std::move(vector); for(auto & mat : matrices) { auto & A = this->getMatrix(mat.first); A.resize(); } std::cout << "Allocating vectors " << solution.get() << std::endl; // should also redo the lumped matrix return ret; } /* -------------------------------------------------------------------------- */ void DOFManagerPETSc::assembleToGlobalArray( const ID & dof_id, const Array & array_to_assemble, SolverVector & global_array, Real scale_factor) { const auto & dof_data = getDOFDataTyped(dof_id); auto & g = dynamic_cast(global_array); AKANTU_DEBUG_ASSERT(array_to_assemble.size() * array_to_assemble.getNbComponent() == dof_data.local_nb_dofs, "The array to assemble does not have the proper size"); g.addValuesLocal(dof_data.local_equation_number_petsc, array_to_assemble, scale_factor); } /* -------------------------------------------------------------------------- */ void DOFManagerPETSc::getArrayPerDOFs(const ID & dof_id, const SolverVector & global_array, Array & local) { const auto & dof_data = getDOFDataTyped(dof_id); const auto & petsc_vector = dynamic_cast(global_array); AKANTU_DEBUG_ASSERT( local.size() * local.getNbComponent() == dof_data.local_nb_dofs, "The array to get the values does not have the proper size"); petsc_vector.getValuesLocal(dof_data.local_equation_number_petsc, local); } /* -------------------------------------------------------------------------- */ void DOFManagerPETSc::assembleElementalMatricesToMatrix( const ID & matrix_id, const ID & dof_id, const Array & elementary_mat, const ElementType & type, const GhostType & ghost_type, const MatrixType & elemental_matrix_type, const Array & filter_elements) { auto & A = getMatrix(matrix_id); DOFManager::assembleElementalMatricesToMatrix_( A, dof_id, elementary_mat, type, ghost_type, elemental_matrix_type, filter_elements); A.applyModifications(); } /* -------------------------------------------------------------------------- */ void DOFManagerPETSc::assemblePreassembledMatrix( const ID & dof_id_m, const ID & dof_id_n, const ID & matrix_id, const TermsToAssemble & terms) { auto & A = getMatrix(matrix_id); DOFManager::assemblePreassembledMatrix_(A, dof_id_m, dof_id_n, terms); A.applyModifications(); } /* -------------------------------------------------------------------------- */ void DOFManagerPETSc::assembleMatMulVectToArray(const ID & dof_id, const ID & A_id, const Array & x, Array & array, Real scale_factor) { DOFManager::assembleMatMulVectToArray_( dof_id, A_id, x, array, scale_factor); } /* -------------------------------------------------------------------------- */ void DOFManagerPETSc::makeConsistentForPeriodicity(const ID & /*dof_id*/, SolverVector & /*array*/) {} /* -------------------------------------------------------------------------- */ NonLinearSolver & DOFManagerPETSc::getNewNonLinearSolver(const ID & id, const NonLinearSolverType & type) { return this->registerNonLinearSolver(*this, id, type); } /* -------------------------------------------------------------------------- */ TimeStepSolver & DOFManagerPETSc::getNewTimeStepSolver(const ID & id, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver) { return this->registerTimeStepSolver(*this, id, type, non_linear_solver); } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManagerPETSc::getNewMatrix(const ID & id, const MatrixType & matrix_type) { return this->registerSparseMatrix(*this, id, matrix_type); } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManagerPETSc::getNewMatrix(const ID & id, const ID & matrix_to_copy_id) { return this->registerSparseMatrix(id, matrix_to_copy_id); } /* -------------------------------------------------------------------------- */ SparseMatrixPETSc & DOFManagerPETSc::getMatrix(const ID & id) { auto & matrix = DOFManager::getMatrix(id); return dynamic_cast(matrix); } /* -------------------------------------------------------------------------- */ SolverVector & DOFManagerPETSc::getNewLumpedMatrix(const ID & id) { return this->registerLumpedMatrix(*this, id); } /* -------------------------------------------------------------------------- */ SolverVectorPETSc & DOFManagerPETSc::getSolution() { return dynamic_cast(*this->solution); } const SolverVectorPETSc & DOFManagerPETSc::getSolution() const { return dynamic_cast(*this->solution); } SolverVectorPETSc & DOFManagerPETSc::getResidual() { return dynamic_cast(*this->residual); } const SolverVectorPETSc & DOFManagerPETSc::getResidual() const { return dynamic_cast(*this->residual); } /* -------------------------------------------------------------------------- */ static bool dof_manager_is_registered [[gnu::unused]] = DOFManagerFactory::getInstance().registerAllocator( "petsc", [](Mesh & mesh, const ID & id, const MemoryID & mem_id) -> std::unique_ptr { return std::make_unique(mesh, id, mem_id); }); } // namespace akantu diff --git a/src/model/common/integration_scheme/integration_scheme_2nd_order.cc b/src/model/common/integration_scheme/integration_scheme_2nd_order.cc index ee1dd2a39..2041da71d 100644 --- a/src/model/common/integration_scheme/integration_scheme_2nd_order.cc +++ b/src/model/common/integration_scheme/integration_scheme_2nd_order.cc @@ -1,106 +1,106 @@ /** * @file integration_scheme_2nd_order.cc * * @author Nicolas Richart * * @date creation: Fri Oct 23 2015 * @date last modification: Wed Jan 31 2018 * * @brief Implementation of the common part of 2nd order integration schemes * * @section LICENSE * * Copyright (©) 2015-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 "integration_scheme_2nd_order.hh" #include "dof_manager.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ std::vector IntegrationScheme2ndOrder::getNeededMatrixList() { return {"K", "M", "C"}; } /* -------------------------------------------------------------------------- */ void IntegrationScheme2ndOrder::predictor(Real delta_t) { AKANTU_DEBUG_IN(); Array & u = this->dof_manager.getDOFs(this->dof_id); Array & u_dot = this->dof_manager.getDOFsDerivatives(this->dof_id, 1); Array & u_dot_dot = this->dof_manager.getDOFsDerivatives(this->dof_id, 2); const Array & blocked_dofs = this->dof_manager.getBlockedDOFs(this->dof_id); this->predictor(delta_t, u, u_dot, u_dot_dot, blocked_dofs); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void IntegrationScheme2ndOrder::corrector(const SolutionType & type, Real delta_t) { AKANTU_DEBUG_IN(); Array & u = this->dof_manager.getDOFs(this->dof_id); Array & u_dot = this->dof_manager.getDOFsDerivatives(this->dof_id, 1); Array & u_dot_dot = this->dof_manager.getDOFsDerivatives(this->dof_id, 2); const Array & solution = this->dof_manager.getSolution(this->dof_id); const Array & blocked_dofs = this->dof_manager.getBlockedDOFs(this->dof_id); this->corrector(type, delta_t, u, u_dot, u_dot_dot, blocked_dofs, solution); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void IntegrationScheme2ndOrder::assembleResidual(bool is_lumped) { AKANTU_DEBUG_IN(); if (this->dof_manager.hasMatrix("C")) { const Array & first_derivative = this->dof_manager.getDOFsDerivatives(this->dof_id, 1); this->dof_manager.assembleMatMulVectToResidual(this->dof_id, "C", first_derivative, -1); } const Array & second_derivative = this->dof_manager.getDOFsDerivatives(this->dof_id, 2); if (not is_lumped) { this->dof_manager.assembleMatMulVectToResidual(this->dof_id, "M", second_derivative, -1); } else { this->dof_manager.assembleLumpedMatMulVectToResidual(this->dof_id, "M", second_derivative, -1); } - + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/common/non_linear_solver/non_linear_solver_petsc.cc b/src/model/common/non_linear_solver/non_linear_solver_petsc.cc index 7b14f5da0..cdeb04ac3 100644 --- a/src/model/common/non_linear_solver/non_linear_solver_petsc.cc +++ b/src/model/common/non_linear_solver/non_linear_solver_petsc.cc @@ -1,206 +1,208 @@ /** * @file non_linear_solver_petsc.cc * * @author Nicolas Richart * * @date creation Mon Dec 31 2018 * * @brief A Documented file. * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "non_linear_solver_petsc.hh" #include "dof_manager_petsc.hh" #include "mpi_communicator_data.hh" #include "solver_callback.hh" #include "solver_vector_petsc.hh" #include "sparse_matrix_petsc.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { NonLinearSolverPETSc::NonLinearSolverPETSc( DOFManagerPETSc & dof_manager, const NonLinearSolverType & non_linear_solver_type, const ID & id, UInt memory_id) : NonLinearSolver(dof_manager, non_linear_solver_type, id, memory_id), dof_manager(dof_manager) { std::unordered_map petsc_non_linear_solver_types{ {NonLinearSolverType::_newton_raphson, SNESNEWTONLS}, {NonLinearSolverType::_linear, SNESKSPONLY}, {NonLinearSolverType::_gmres, SNESNGMRES}, {NonLinearSolverType::_bfgs, SNESQN}, {NonLinearSolverType::_cg, SNESNCG}}; this->has_internal_set_param = true; for (const auto & pair : petsc_non_linear_solver_types) { supported_type.insert(pair.first); } this->checkIfTypeIsSupported(); auto mpi_comm = dof_manager.getMPIComm(); PETSc_call(SNESCreate, mpi_comm, &snes); auto it = petsc_non_linear_solver_types.find(non_linear_solver_type); if (it != petsc_non_linear_solver_types.end()) { PETSc_call(SNESSetType, snes, it->second); } SNESSetFromOptions(snes); } /* -------------------------------------------------------------------------- */ NonLinearSolverPETSc::~NonLinearSolverPETSc() { PETSc_call(SNESDestroy, &snes); } /* -------------------------------------------------------------------------- */ class NonLinearSolverPETScCallback { public: NonLinearSolverPETScCallback(DOFManagerPETSc & dof_manager, SolverVectorPETSc & x) : dof_manager(dof_manager), x(x), x_prev(x, "previous_solution") {} void corrector() { auto & dx = dof_manager.getSolution(); PETSc_call(VecWAXPY, dx, -1., x_prev, x); - VecView(x_prev, PETSC_VIEWER_STDOUT_WORLD); - VecView(x, PETSC_VIEWER_STDOUT_WORLD); dof_manager.splitSolutionPerDOFs(); callback->corrector(); PETSc_call(VecCopy, x, x_prev); } void assembleResidual() { corrector(); callback->assembleResidual(); - - auto & r = - dynamic_cast(dof_manager.getResidual()); - VecView(r, PETSC_VIEWER_STDOUT_WORLD); } void assembleJacobian() { //corrector(); callback->assembleMatrix("J"); } void setInitialSolution(SolverVectorPETSc & x) { PETSc_call(VecCopy, x, x_prev); } void setCallback(SolverCallback & callback) { this->callback = &callback; } private: // SNES & snes; SolverCallback * callback; DOFManagerPETSc & dof_manager; SolverVectorPETSc & x; SolverVectorPETSc x_prev; }; // namespace akantu /* -------------------------------------------------------------------------- */ PetscErrorCode NonLinearSolverPETSc::FormFunction(SNES /*snes*/, Vec /*dx*/, Vec /*f*/, void * ctx) { auto * _this = reinterpret_cast(ctx); _this->assembleResidual(); return 0; } /* -------------------------------------------------------------------------- */ PetscErrorCode NonLinearSolverPETSc::FormJacobian(SNES /*snes*/, Vec /*dx*/, Mat /*J*/, Mat /*P*/, void * ctx) { auto * _this = reinterpret_cast(ctx); _this->assembleJacobian(); return 0; } /* -------------------------------------------------------------------------- */ void NonLinearSolverPETSc::solve(SolverCallback & callback) { this->dof_manager.updateGlobalBlockedDofs(); callback.assembleMatrix("J"); - auto & x_ = dof_manager.getSolution(); - - if (not x or x->size() != x_.size()) { - x = std::make_unique(x_, "temporary_solution"); + auto & global_x = dof_manager.getSolution(); + global_x.clear(); + + if (not x) { + x = std::make_unique(global_x, "temporary_solution"); } - x->clear(); + *x = global_x; + if (not ctx) { ctx = std::make_unique(dof_manager, *x); } + ctx->setCallback(callback); - ctx->setInitialSolution(x_); - + ctx->setInitialSolution(global_x); + auto & rhs = dof_manager.getResidual(); auto & J = dof_manager.getMatrix("J"); + PETSc_call(SNESSetFunction, snes, rhs, NonLinearSolverPETSc::FormFunction, ctx.get()); PETSc_call(SNESSetJacobian, snes, J, J, NonLinearSolverPETSc::FormJacobian, ctx.get()); rhs.clear(); callback.predictor(); + callback.assembleResidual(); + PETSc_call(SNESSolve, snes, nullptr, *x); - PETSc_call(VecAXPY, x_, -1.0, *x); + PETSc_call(VecAXPY, global_x, -1.0, *x); + + dof_manager.splitSolutionPerDOFs(); callback.corrector(); } /* -------------------------------------------------------------------------- */ void NonLinearSolverPETSc::set_param(const ID & param, const std::string & value) { std::map akantu_to_petsc_option = {{"max_iterations", "snes_max_it"}, {"threshold", "snes_stol"}}; auto it = akantu_to_petsc_option.find(param); auto p = it == akantu_to_petsc_option.end() ? param : it->second; PetscOptionsSetValue(NULL, p.c_str(), value.c_str()); SNESSetFromOptions(snes); PetscOptionsClear(NULL); } /* -------------------------------------------------------------------------- */ void NonLinearSolverPETSc::parseSection(const ParserSection & section) { auto parameters = section.getParameters(); for (auto && param : range(parameters.first, parameters.second)) { PetscOptionsSetValue(NULL, param.getName().c_str(), param.getValue().c_str()); } SNESSetFromOptions(snes); PetscOptionsClear(NULL); } } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc index 32e9e8000..260a914d2 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model.cc @@ -1,1195 +1,1195 @@ /** * @file solid_mechanics_model.cc * * @author Ramin Aghababaei * @author Guillaume Anciaux * @author Aurelia Isabel Cuba Ramos * @author David Simon Kammer * @author Daniel Pino Muñoz * @author Nicolas Richart * @author Clement Roux * @author Marco Vocialta * * @date creation: Tue Jul 27 2010 * @date last modification: Wed Feb 21 2018 * * @brief Implementation of the SolidMechanicsModel 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 "solid_mechanics_model.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #include "solid_mechanics_model_tmpl.hh" #include "communicator.hh" #include "element_synchronizer.hh" #include "sparse_matrix.hh" #include "synchronizer_registry.hh" #include "dumpable_inline_impl.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_iohelper_paraview.hh" #endif #include "material_non_local.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /** * A solid mechanics model need a mesh and a dimension to be created. the model * by it self can not do a lot, the good init functions should be called in * order to configure the model depending on what we want to do. * * @param mesh mesh representing the model we want to simulate * @param dim spatial dimension of the problem, if dim = 0 (default value) the * dimension of the problem is assumed to be the on of the mesh * @param id an id to identify the model */ SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id, const ModelType model_type) : Model(mesh, model_type, dim, id, memory_id), BoundaryCondition(), f_m2a(1.0), displacement(nullptr), previous_displacement(nullptr), displacement_increment(nullptr), mass(nullptr), velocity(nullptr), acceleration(nullptr), external_force(nullptr), internal_force(nullptr), blocked_dofs(nullptr), current_position(nullptr), material_index("material index", id, memory_id), material_local_numbering("material local numbering", id, memory_id), - increment_flag(false), are_materials_instantiated(false) { + are_materials_instantiated(false) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("SolidMechanicsFEEngine", mesh, Model::spatial_dimension); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("paraview_all", id, true); this->mesh.addDumpMesh(mesh, Model::spatial_dimension, _not_ghost, _ek_regular); #endif material_selector = std::make_shared(material_index), this->initDOFManager(); this->registerDataAccessor(*this); if (this->mesh.isDistributed()) { auto & synchronizer = this->mesh.getElementSynchronizer(); this->registerSynchronizer(synchronizer, SynchronizationTag::_material_id); this->registerSynchronizer(synchronizer, SynchronizationTag::_smm_mass); this->registerSynchronizer(synchronizer, SynchronizationTag::_smm_stress); this->registerSynchronizer(synchronizer, SynchronizationTag::_for_dump); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModel::~SolidMechanicsModel() { AKANTU_DEBUG_IN(); for (auto & internal : this->registered_internals) { delete internal.second; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setTimeStep(Real time_step, const ID & solver_id) { Model::setTimeStep(time_step, solver_id); #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper().setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ /* Initialization */ /* -------------------------------------------------------------------------- */ /** * This function groups many of the initialization in on function. For most of * basics case the function should be enough. The functions initialize the * model, the internal vectors, set them to 0, and depending on the parameters * it also initialize the explicit or implicit solver. * * @param material_file the file containing the materials to use * @param method the analysis method wanted. See the akantu::AnalysisMethod for * the different possibilities */ void SolidMechanicsModel::initFullImpl(const ModelOptions & options) { material_index.initialize(mesh, _element_kind = _ek_not_defined, _default_value = UInt(-1), _with_nb_element = true); material_local_numbering.initialize(mesh, _element_kind = _ek_not_defined, _with_nb_element = true); Model::initFullImpl(options); // initialize the materials if (this->parser.getLastParsedFile() != "") { this->instantiateMaterials(); } this->initMaterials(); this->initBC(*this, *displacement, *displacement_increment, *external_force); } /* -------------------------------------------------------------------------- */ TimeStepSolverType SolidMechanicsModel::getDefaultSolverType() const { return TimeStepSolverType::_dynamic_lumped; } /* -------------------------------------------------------------------------- */ ModelSolverOptions SolidMechanicsModel::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::_static: { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } case TimeStepSolverType::_dynamic: { if (this->method == _explicit_consistent_mass) { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; } else { options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_trapezoidal_rule_2; options.solution_type["displacement"] = IntegrationScheme::_displacement; } break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); } return options; } /* -------------------------------------------------------------------------- */ std::tuple SolidMechanicsModel::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); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType) { auto & dof_manager = this->getDOFManager(); /* ------------------------------------------------------------------------ */ // for alloc type of solvers this->allocNodalField(this->displacement, spatial_dimension, "displacement"); this->allocNodalField(this->previous_displacement, spatial_dimension, "previous_displacement"); this->allocNodalField(this->displacement_increment, spatial_dimension, "displacement_increment"); this->allocNodalField(this->internal_force, spatial_dimension, "internal_force"); this->allocNodalField(this->external_force, spatial_dimension, "external_force"); this->allocNodalField(this->blocked_dofs, spatial_dimension, "blocked_dofs"); this->allocNodalField(this->current_position, spatial_dimension, "current_position"); // initialize the current positions this->current_position->copy(this->mesh.getNodes()); /* ------------------------------------------------------------------------ */ if (!dof_manager.hasDOFs("displacement")) { dof_manager.registerDOFs("displacement", *this->displacement, _dst_nodal); dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs); dof_manager.registerDOFsIncrement("displacement", *this->displacement_increment); dof_manager.registerDOFsPrevious("displacement", *this->previous_displacement); } /* ------------------------------------------------------------------------ */ // for dynamic if (time_step_solver_type == TimeStepSolverType::_dynamic || time_step_solver_type == TimeStepSolverType::_dynamic_lumped) { this->allocNodalField(this->velocity, spatial_dimension, "velocity"); this->allocNodalField(this->acceleration, spatial_dimension, "acceleration"); if (!dof_manager.hasDOFsDerivatives("displacement", 1)) { dof_manager.registerDOFsDerivative("displacement", 1, *this->velocity); dof_manager.registerDOFsDerivative("displacement", 2, *this->acceleration); } } } /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian */ void SolidMechanicsModel::initModel() { /// \todo add the current position as a parameter to initShapeFunctions for /// large deformation getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleResidual() { AKANTU_DEBUG_IN(); /* ------------------------------------------------------------------------ */ // computes the internal forces this->assembleInternalForces(); /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", *this->external_force, 1); this->getDOFManager().assembleToResidual("displacement", *this->internal_force, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleResidual(const ID & residual_part) { AKANTU_DEBUG_IN(); if ("external" == residual_part) { this->getDOFManager().assembleToResidual("displacement", *this->external_force, 1); AKANTU_DEBUG_OUT(); return; } if ("internal" == residual_part) { this->getDOFManager().assembleToResidual("displacement", *this->internal_force, 1); AKANTU_DEBUG_OUT(); return; } AKANTU_CUSTOM_EXCEPTION( debug::SolverCallbackResidualPartUnknown(residual_part)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MatrixType SolidMechanicsModel::getMatrixType(const ID & matrix_id) { // \TODO check the materials to know what is the correct answer if (matrix_id == "C") return _mt_not_defined; return _symmetric; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { this->assembleStiffnessMatrix(); } else if (matrix_id == "M") { this->assembleMass(); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleLumpedMatrix(const ID & matrix_id) { if (matrix_id == "M") { this->assembleMassLumped(); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::beforeSolveStep() { for (auto & material : materials) material->beforeSolveStep(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::afterSolveStep() { for (auto & material : materials) material->afterSolveStep(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::predictor() { ++displacement_release; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::corrector() { ++displacement_release; } /* -------------------------------------------------------------------------- */ /** * This function computes the internal forces as F_{int} = \int_{\Omega} N * \sigma d\Omega@f$ */ void SolidMechanicsModel::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); this->internal_force->clear(); // compute the stresses of local elements AKANTU_DEBUG_INFO("Compute local stresses"); for (auto & material : materials) { material->computeAllStresses(_not_ghost); } /* ------------------------------------------------------------------------ */ /* Computation of the non local part */ if (this->non_local_manager) this->non_local_manager->computeAllNonLocalStresses(); // communicate the stresses AKANTU_DEBUG_INFO("Send data for residual assembly"); this->asynchronousSynchronize(SynchronizationTag::_smm_stress); // assemble the forces due to local stresses AKANTU_DEBUG_INFO("Assemble residual for local elements"); for (auto & material : materials) { material->assembleInternalForces(_not_ghost); } // finalize communications AKANTU_DEBUG_INFO("Wait distant stresses"); this->waitEndSynchronize(SynchronizationTag::_smm_stress); // assemble the stresses due to ghost elements AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); for (auto & material : materials) { material->assembleInternalForces(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix."); // Check if materials need to recompute the matrix bool need_to_reassemble = false; for (auto & material : materials) { need_to_reassemble |= material->hasStiffnessMatrixChanged(); } if (need_to_reassemble) { this->getDOFManager().getMatrix("K").clear(); // call compute stiffness matrix on each local elements for (auto & material : materials) { material->assembleStiffnessMatrix(_not_ghost); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateCurrentPosition() { if (this->current_position_release == this->displacement_release) return; this->current_position->copy(this->mesh.getNodes()); auto cpos_it = this->current_position->begin(Model::spatial_dimension); auto cpos_end = this->current_position->end(Model::spatial_dimension); auto disp_it = this->displacement->begin(Model::spatial_dimension); for (; cpos_it != cpos_end; ++cpos_it, ++disp_it) { *cpos_it += *disp_it; } this->current_position_release = this->displacement_release; } /* -------------------------------------------------------------------------- */ const Array & SolidMechanicsModel::getCurrentPosition() { this->updateCurrentPosition(); return *this->current_position; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateDataForNonLocalCriterion( ElementTypeMapReal & criterion) { const ID field_name = criterion.getName(); for (auto & material : materials) { if (!material->isInternal(field_name, _ek_regular)) continue; for (auto ghost_type : ghost_types) { material->flattenInternal(field_name, criterion, ghost_type, _ek_regular); } } } /* -------------------------------------------------------------------------- */ /* Information */ /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep() { AKANTU_DEBUG_IN(); Real min_dt = getStableTimeStep(_not_ghost); /// reduction min over all processors mesh.getCommunicator().allReduce(min_dt, SynchronizerOperation::_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); Real min_dt = std::numeric_limits::max(); this->updateCurrentPosition(); Element elem; elem.ghost_type = ghost_type; for (auto type : mesh.elementTypes(Model::spatial_dimension, ghost_type, _ek_regular)) { elem.type = type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); UInt nb_element = mesh.getNbElement(type); auto mat_indexes = material_index(type, ghost_type).begin(); auto mat_loc_num = material_local_numbering(type, ghost_type).begin(); Array X(0, nb_nodes_per_element * Model::spatial_dimension); FEEngine::extractNodalToElementField(mesh, *current_position, X, type, _not_ghost); auto X_el = X.begin(Model::spatial_dimension, nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++X_el, ++mat_indexes, ++mat_loc_num) { elem.element = *mat_loc_num; Real el_h = getFEEngine().getElementInradius(*X_el, type); Real el_c = this->materials[*mat_indexes]->getCelerity(elem); Real el_dt = el_h / el_c; min_dt = std::min(min_dt, el_dt); } } AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy() { AKANTU_DEBUG_IN(); Real ekin = 0.; UInt nb_nodes = mesh.getNbNodes(); if (this->getDOFManager().hasLumpedMatrix("M")) { auto m_it = this->mass->begin(Model::spatial_dimension); auto m_end = this->mass->end(Model::spatial_dimension); auto v_it = this->velocity->begin(Model::spatial_dimension); for (UInt n = 0; m_it != m_end; ++n, ++m_it, ++v_it) { const auto & v = *v_it; const auto & m = *m_it; Real mv2 = 0.; auto is_local_node = mesh.isLocalOrMasterNode(n); // bool is_not_pbc_slave_node = !isPBCSlaveNode(n); auto count_node = is_local_node; // && is_not_pbc_slave_node; if (count_node) { for (UInt i = 0; i < Model::spatial_dimension; ++i) { if (m(i) > std::numeric_limits::epsilon()) mv2 += v(i) * v(i) * m(i); } } ekin += mv2; } } else if (this->getDOFManager().hasMatrix("M")) { Array Mv(nb_nodes, Model::spatial_dimension); this->getDOFManager().assembleMatMulVectToArray("displacement", "M", *this->velocity, Mv); for (auto && data : zip(arange(nb_nodes), make_view(Mv, spatial_dimension), make_view(*this->velocity, spatial_dimension))) { ekin += std::get<2>(data).dot(std::get<1>(data)) * mesh.isLocalOrMasterNode(std::get<0>(data)); } } else { AKANTU_ERROR("No function called to assemble the mass matrix."); } mesh.getCommunicator().allReduce(ekin, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return ekin * .5; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy(const ElementType & type, UInt index) { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type); Array vel_on_quad(nb_quadrature_points, Model::spatial_dimension); Array filter_element(1, 1, index); getFEEngine().interpolateOnIntegrationPoints(*velocity, vel_on_quad, Model::spatial_dimension, type, _not_ghost, filter_element); auto vit = vel_on_quad.begin(Model::spatial_dimension); auto vend = vel_on_quad.end(Model::spatial_dimension); Vector rho_v2(nb_quadrature_points); Real rho = materials[material_index(type)(index)]->getRho(); for (UInt q = 0; vit != vend; ++vit, ++q) { rho_v2(q) = rho * vit->dot(*vit); } AKANTU_DEBUG_OUT(); return .5 * getFEEngine().integrate(rho_v2, type, index); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getExternalWork() { AKANTU_DEBUG_IN(); auto ext_force_it = external_force->begin(Model::spatial_dimension); auto int_force_it = internal_force->begin(Model::spatial_dimension); auto boun_it = blocked_dofs->begin(Model::spatial_dimension); decltype(ext_force_it) incr_or_velo_it; if (this->method == _static) { incr_or_velo_it = this->displacement_increment->begin(Model::spatial_dimension); } else { incr_or_velo_it = this->velocity->begin(Model::spatial_dimension); } Real work = 0.; UInt nb_nodes = this->mesh.getNbNodes(); for (UInt n = 0; n < nb_nodes; ++n, ++ext_force_it, ++int_force_it, ++boun_it, ++incr_or_velo_it) { const auto & int_force = *int_force_it; const auto & ext_force = *ext_force_it; const auto & boun = *boun_it; const auto & incr_or_velo = *incr_or_velo_it; bool is_local_node = this->mesh.isLocalOrMasterNode(n); // bool is_not_pbc_slave_node = !this->isPBCSlaveNode(n); bool count_node = is_local_node; // && is_not_pbc_slave_node; if (count_node) { for (UInt i = 0; i < Model::spatial_dimension; ++i) { if (boun(i)) work -= int_force(i) * incr_or_velo(i); else work += ext_force(i) * incr_or_velo(i); } } } mesh.getCommunicator().allReduce(work, SynchronizerOperation::_sum); if (this->method != _static) work *= this->getTimeStep(); AKANTU_DEBUG_OUT(); return work; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id) { AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(); } else if (energy_id == "external work") { return getExternalWork(); } Real energy = 0.; for (auto & material : materials) energy += material->getEnergy(energy_id); /// reduction sum over all processors mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id, const ElementType & type, UInt index) { AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(type, index); } UInt mat_index = this->material_index(type, _not_ghost)(index); UInt mat_loc_num = this->material_local_numbering(type, _not_ghost)(index); Real energy = this->materials[mat_index]->getEnergy(energy_id, type, mat_loc_num); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsAdded(const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); this->material_index.initialize(mesh, _element_kind = _ek_not_defined, _with_nb_element = true, _default_value = UInt(-1)); this->material_local_numbering.initialize( mesh, _element_kind = _ek_not_defined, _with_nb_element = true, _default_value = UInt(-1)); ElementTypeMapArray filter("new_element_filter", this->getID(), this->getMemoryID()); for (auto & elem : element_list) { if (!filter.exists(elem.type, elem.ghost_type)) filter.alloc(0, 1, elem.type, elem.ghost_type); filter(elem.type, elem.ghost_type).push_back(elem.element); } this->assignMaterialToElements(&filter); for (auto & material : materials) material->onElementsAdded(element_list, event); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsRemoved( const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) { for (auto & material : materials) { material->onElementsRemoved(element_list, new_numbering, event); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); if (displacement) { displacement->resize(nb_nodes, 0.); ++displacement_release; } if (mass) mass->resize(nb_nodes, 0.); if (velocity) velocity->resize(nb_nodes, 0.); if (acceleration) acceleration->resize(nb_nodes, 0.); if (external_force) external_force->resize(nb_nodes, 0.); if (internal_force) internal_force->resize(nb_nodes, 0.); if (blocked_dofs) blocked_dofs->resize(nb_nodes, 0.); if (current_position) current_position->resize(nb_nodes, 0.); if (previous_displacement) previous_displacement->resize(nb_nodes, 0.); if (displacement_increment) displacement_increment->resize(nb_nodes, 0.); for (auto & material : materials) { material->onNodesAdded(nodes_list, event); } need_to_reassemble_lumped_mass = true; need_to_reassemble_mass = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesRemoved(const Array & /*element_list*/, const Array & new_numbering, const RemovedNodesEvent & /*event*/) { if (displacement) { mesh.removeNodesFromArray(*displacement, new_numbering); ++displacement_release; } if (mass) mesh.removeNodesFromArray(*mass, new_numbering); if (velocity) mesh.removeNodesFromArray(*velocity, new_numbering); if (acceleration) mesh.removeNodesFromArray(*acceleration, new_numbering); if (internal_force) mesh.removeNodesFromArray(*internal_force, new_numbering); if (external_force) mesh.removeNodesFromArray(*external_force, new_numbering); if (blocked_dofs) mesh.removeNodesFromArray(*blocked_dofs, new_numbering); // if (increment_acceleration) // mesh.removeNodesFromArray(*increment_acceleration, new_numbering); if (displacement_increment) mesh.removeNodesFromArray(*displacement_increment, new_numbering); if (previous_displacement) mesh.removeNodesFromArray(*previous_displacement, new_numbering); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space(indent, AKANTU_INDENT); stream << space << "Solid Mechanics 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 << " ]" << std::endl; stream << space << " + nodals information [" << std::endl; displacement->printself(stream, indent + 2); if (velocity) velocity->printself(stream, indent + 2); if (acceleration) acceleration->printself(stream, indent + 2); if (mass) mass->printself(stream, indent + 2); external_force->printself(stream, indent + 2); internal_force->printself(stream, indent + 2); blocked_dofs->printself(stream, indent + 2); stream << space << " ]" << std::endl; stream << space << " + material information [" << std::endl; material_index.printself(stream, indent + 2); stream << space << " ]" << std::endl; stream << space << " + materials [" << std::endl; for (auto & material : materials) material->printself(stream, indent + 2); stream << space << " ]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initializeNonLocal() { this->non_local_manager->synchronize(*this, SynchronizationTag::_material_id); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::insertIntegrationPointsInNeighborhoods( const GhostType & ghost_type) { for (auto & mat : materials) { MaterialNonLocalInterface * mat_non_local; if ((mat_non_local = dynamic_cast(mat.get())) == nullptr) continue; ElementTypeMapArray quadrature_points_coordinates( "quadrature_points_coordinates_tmp_nl", this->id, this->memory_id); quadrature_points_coordinates.initialize(this->getFEEngine(), _nb_component = spatial_dimension, _ghost_type = ghost_type); for (auto & type : quadrature_points_coordinates.elementTypes( Model::spatial_dimension, ghost_type)) { this->getFEEngine().computeIntegrationPointsCoordinates( quadrature_points_coordinates(type, ghost_type), type, ghost_type); } mat_non_local->initMaterialNonLocal(); mat_non_local->insertIntegrationPointsInNeighborhoods( ghost_type, quadrature_points_coordinates); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeNonLocalStresses( const GhostType & ghost_type) { for (auto & mat : materials) { if (dynamic_cast(mat.get()) == nullptr) continue; auto & mat_non_local = dynamic_cast(*mat); mat_non_local.computeNonLocalStresses(ghost_type); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateLocalInternal( ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind) { const ID field_name = internal_flat.getName(); for (auto & material : materials) { if (material->isInternal(field_name, kind)) material->flattenInternal(field_name, internal_flat, ghost_type, kind); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateNonLocalInternal( ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind) { const ID field_name = internal_flat.getName(); for (auto & mat : materials) { if (dynamic_cast(mat.get()) == nullptr) continue; auto & mat_non_local = dynamic_cast(*mat); mat_non_local.updateNonLocalInternals(internal_flat, field_name, ghost_type, kind); } } /* -------------------------------------------------------------------------- */ FEEngine & SolidMechanicsModel::getFEEngineBoundary(const ID & name) { return dynamic_cast( getFEEngineClassBoundary(name)); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::splitElementByMaterial( const Array & elements, std::vector> & elements_per_mat) const { for (const auto & el : elements) { Element mat_el = el; mat_el.element = this->material_local_numbering(el); elements_per_mat[this->material_index(el)].push_back(mat_el); } } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModel::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::_material_id: { size += elements.size() * sizeof(UInt); break; } case SynchronizationTag::_smm_mass: { size += nb_nodes_per_element * sizeof(Real) * Model::spatial_dimension; // mass vector break; } case SynchronizationTag::_smm_for_gradu: { size += nb_nodes_per_element * Model::spatial_dimension * sizeof(Real); // displacement break; } case SynchronizationTag::_smm_boundary: { // force, displacement, boundary size += nb_nodes_per_element * Model::spatial_dimension * (2 * sizeof(Real) + sizeof(bool)); break; } case SynchronizationTag::_for_dump: { // displacement, velocity, acceleration, residual, force size += nb_nodes_per_element * Model::spatial_dimension * sizeof(Real) * 5; break; } default: {} } if (tag != SynchronizationTag::_material_id) { splitByMaterial(elements, [&](auto && mat, auto && elements) { size += mat.getNbData(elements, tag); }); } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); switch (tag) { case SynchronizationTag::_material_id: { this->packElementalDataHelper(material_index, buffer, elements, false, getFEEngine()); break; } case SynchronizationTag::_smm_mass: { packNodalDataHelper(*mass, buffer, elements, mesh); break; } case SynchronizationTag::_smm_for_gradu: { packNodalDataHelper(*displacement, buffer, elements, mesh); break; } case SynchronizationTag::_for_dump: { packNodalDataHelper(*displacement, buffer, elements, mesh); packNodalDataHelper(*velocity, buffer, elements, mesh); packNodalDataHelper(*acceleration, buffer, elements, mesh); packNodalDataHelper(*internal_force, buffer, elements, mesh); packNodalDataHelper(*external_force, buffer, elements, mesh); break; } case SynchronizationTag::_smm_boundary: { packNodalDataHelper(*external_force, buffer, elements, mesh); packNodalDataHelper(*velocity, buffer, elements, mesh); packNodalDataHelper(*blocked_dofs, buffer, elements, mesh); break; } default: {} } if (tag != SynchronizationTag::_material_id) { splitByMaterial(elements, [&](auto && mat, auto && elements) { mat.packData(buffer, elements, tag); }); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); switch (tag) { case SynchronizationTag::_material_id: { for (auto && element : elements) { UInt recv_mat_index; buffer >> recv_mat_index; UInt & mat_index = material_index(element); if (mat_index != UInt(-1)) continue; // add ghosts element to the correct material mat_index = recv_mat_index; UInt index = materials[mat_index]->addElement(element); material_local_numbering(element) = index; } break; } case SynchronizationTag::_smm_mass: { unpackNodalDataHelper(*mass, buffer, elements, mesh); break; } case SynchronizationTag::_smm_for_gradu: { unpackNodalDataHelper(*displacement, buffer, elements, mesh); break; } case SynchronizationTag::_for_dump: { unpackNodalDataHelper(*displacement, buffer, elements, mesh); unpackNodalDataHelper(*velocity, buffer, elements, mesh); unpackNodalDataHelper(*acceleration, buffer, elements, mesh); unpackNodalDataHelper(*internal_force, buffer, elements, mesh); unpackNodalDataHelper(*external_force, buffer, elements, mesh); break; } case SynchronizationTag::_smm_boundary: { unpackNodalDataHelper(*external_force, buffer, elements, mesh); unpackNodalDataHelper(*velocity, buffer, elements, mesh); unpackNodalDataHelper(*blocked_dofs, buffer, elements, mesh); break; } default: {} } if (tag != SynchronizationTag::_material_id) { splitByMaterial(elements, [&](auto && mat, auto && elements) { mat.unpackData(buffer, elements, tag); }); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModel::getNbData(const Array & dofs, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; // UInt nb_nodes = mesh.getNbNodes(); switch (tag) { case SynchronizationTag::_smm_uv: { size += sizeof(Real) * Model::spatial_dimension * 2; break; } case SynchronizationTag::_smm_res: { size += sizeof(Real) * Model::spatial_dimension; break; } case SynchronizationTag::_smm_mass: { size += sizeof(Real) * Model::spatial_dimension; break; } case SynchronizationTag::_for_dump: { size += sizeof(Real) * Model::spatial_dimension * 5; break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size * dofs.size(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); switch (tag) { case SynchronizationTag::_smm_uv: { packDOFDataHelper(*displacement, buffer, dofs); packDOFDataHelper(*velocity, buffer, dofs); break; } case SynchronizationTag::_smm_res: { packDOFDataHelper(*internal_force, buffer, dofs); break; } case SynchronizationTag::_smm_mass: { packDOFDataHelper(*mass, buffer, dofs); break; } case SynchronizationTag::_for_dump: { packDOFDataHelper(*displacement, buffer, dofs); packDOFDataHelper(*velocity, buffer, dofs); packDOFDataHelper(*acceleration, buffer, dofs); packDOFDataHelper(*internal_force, buffer, dofs); packDOFDataHelper(*external_force, buffer, dofs); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); switch (tag) { case SynchronizationTag::_smm_uv: { unpackDOFDataHelper(*displacement, buffer, dofs); unpackDOFDataHelper(*velocity, buffer, dofs); break; } case SynchronizationTag::_smm_res: { unpackDOFDataHelper(*internal_force, buffer, dofs); break; } case SynchronizationTag::_smm_mass: { unpackDOFDataHelper(*mass, buffer, dofs); break; } case SynchronizationTag::_for_dump: { unpackDOFDataHelper(*displacement, buffer, dofs); unpackDOFDataHelper(*velocity, buffer, dofs); unpackDOFDataHelper(*acceleration, buffer, dofs); unpackDOFDataHelper(*internal_force, buffer, dofs); unpackDOFDataHelper(*external_force, buffer, dofs); break; } default: { AKANTU_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh index 0edab13c7..a568e4a8f 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,568 +1,565 @@ /** * @file solid_mechanics_model.hh * * @author Guillaume Anciaux * @author Daniel Pino Muñoz * @author Nicolas Richart * * @date creation: Tue Jul 27 2010 * @date last modification: Wed Feb 21 2018 * * @brief Model of Solid Mechanics * * @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 "model.hh" #include "non_local_manager.hh" #include "solid_mechanics_model_event_handler.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_MECHANICS_MODEL_HH__ #define __AKANTU_SOLID_MECHANICS_MODEL_HH__ namespace akantu { class Material; class MaterialSelector; class DumperIOHelper; class NonLocalManager; template class IntegratorGauss; template class ShapeLagrange; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ class SolidMechanicsModel : public Model, public DataAccessor, public DataAccessor, public BoundaryCondition, public NonLocalManagerCallback, public EventHandlerManager { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: class NewMaterialElementsEvent : public NewElementsEvent { public: AKANTU_GET_MACRO_NOT_CONST(MaterialList, material, Array &); AKANTU_GET_MACRO(MaterialList, material, const Array &); protected: Array material; }; using MyFEEngineType = FEEngineTemplate; protected: using EventManager = EventHandlerManager; public: SolidMechanicsModel( Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "solid_mechanics_model", const MemoryID & memory_id = 0, const ModelType model_type = ModelType::_solid_mechanics_model); ~SolidMechanicsModel() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ protected: /// initialize completely the model void initFullImpl( const ModelOptions & options = SolidMechanicsModelOptions()) override; /// initialize all internal arrays for materials virtual void initMaterials(); /// initialize the model void initModel() override; /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; /* ------------------------------------------------------------------------ */ /* Solver interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the stiffness matrix, virtual void assembleStiffnessMatrix(); /// assembles the internal forces in the array internal_forces virtual void assembleInternalForces(); protected: /// callback for the solver, this adds f_{ext} - f_{int} to the residual void assembleResidual() override; /// callback for the solver, this adds f_{ext} or f_{int} to the residual void assembleResidual(const ID & residual_part) override; bool canSplitResidual() 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 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() override; /// Callback for the model to instantiate the matricees when needed void initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) override; protected: /* ------------------------------------------------------------------------ */ TimeStepSolverType getDefaultSolverType() const override; /* ------------------------------------------------------------------------ */ ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; public: bool isDefaultSolverExplicit() { return method == _explicit_lumped_mass || method == _explicit_consistent_mass; } protected: /// update the current position vector void updateCurrentPosition(); /* ------------------------------------------------------------------------ */ /* Materials (solid_mechanics_model_material.cc) */ /* ------------------------------------------------------------------------ */ public: /// register an empty material of a given type Material & registerNewMaterial(const ID & mat_name, const ID & mat_type, const ID & opt_param); /// reassigns materials depending on the material selector virtual void reassignMaterial(); /// apply a constant eigen_grad_u on all quadrature points of a given material virtual void applyEigenGradU(const Matrix & prescribed_eigen_grad_u, const ID & material_name, const GhostType ghost_type = _not_ghost); protected: /// register a material in the dynamic database Material & registerNewMaterial(const ParserSection & mat_section); /// read the material files to instantiate all the materials void instantiateMaterials(); /// set the element_id_by_material and add the elements to the good materials virtual void assignMaterialToElements(const ElementTypeMapArray * filter = nullptr); /* ------------------------------------------------------------------------ */ /* Mass (solid_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the mass matrix for consistent mass resolutions void assembleMass(); 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); /// fill a vector of rho void computeRho(Array & rho, ElementType type, GhostType ghost_type); /// compute the kinetic energy Real getKineticEnergy(); Real getKineticEnergy(const ElementType & type, UInt index); /// compute the external work (for impose displacement, the velocity should be /// given too) Real getExternalWork(); /* ------------------------------------------------------------------------ */ /* NonLocalManager inherited members */ /* ------------------------------------------------------------------------ */ protected: void initializeNonLocal() override; void updateDataForNonLocalCriterion(ElementTypeMapReal & criterion) override; void computeNonLocalStresses(const GhostType & ghost_type) override; void insertIntegrationPointsInNeighborhoods(const GhostType & ghost_type) override; /// update the values of the non local internal void updateLocalInternal(ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind) override; /// copy the results of the averaging in the materials void updateNonLocalInternal(ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind) override; /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; UInt getNbData(const Array & dofs, const SynchronizationTag & tag) const override; void packData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) const override; void unpackData(CommunicationBuffer & buffer, const Array & dofs, const SynchronizationTag & tag) override; protected: void splitElementByMaterial(const Array & elements, std::vector> & elements_per_mat) const; template void splitByMaterial(const Array & elements, Operation && op) const; /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) override; void onNodesRemoved(const Array & element_list, const Array & new_numbering, const RemovedNodesEvent & event) override; void onElementsAdded(const Array & nodes_list, const NewElementsEvent & event) override; void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) override; void onElementsChanged(const Array &, const Array &, const ElementTypeMapArray &, const ChangedElementsEvent &) override{}; /* ------------------------------------------------------------------------ */ /* Dumpable interface (kept for convenience) and dumper relative functions */ /* ------------------------------------------------------------------------ */ public: virtual void onDump(); //! decide wether a field is a material internal or not bool isInternal(const std::string & field_name, const ElementKind & element_kind); #ifndef SWIG //! give the amount of data per element virtual ElementTypeMap getInternalDataPerElem(const std::string & field_name, const ElementKind & kind); //! flatten a given material internal field ElementTypeMapArray & flattenInternal(const std::string & field_name, const ElementKind & kind, const GhostType ghost_type = _not_ghost); //! flatten all the registered material internals void flattenAllRegisteredInternals(const ElementKind & kind); #endif dumper::Field * createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; dumper::Field * createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; dumper::Field * createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const 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); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt); /// set the value of the time step void setTimeStep(Real time_step, const ID & solver_id = "") override; /// get the value of the conversion from forces/ mass to acceleration AKANTU_GET_MACRO(F_M2A, f_m2a, Real); /// set the value of the conversion from forces/ mass to acceleration AKANTU_SET_MACRO(F_M2A, f_m2a, Real); /// get the SolidMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement, Array &); /// get the SolidMechanicsModel::previous_displacement vector AKANTU_GET_MACRO(PreviousDisplacement, *previous_displacement, Array &); /// get the SolidMechanicsModel::current_position vector \warn only consistent /// after a call to SolidMechanicsModel::updateCurrentPosition const Array & getCurrentPosition(); /// get the SolidMechanicsModel::increment vector \warn only consistent if /// SolidMechanicsModel::setIncrementFlagOn has been called before AKANTU_GET_MACRO(Increment, *displacement_increment, Array &); /// get the lumped SolidMechanicsModel::mass vector AKANTU_GET_MACRO(Mass, *mass, Array &); /// get the SolidMechanicsModel::velocity vector AKANTU_GET_MACRO(Velocity, *velocity, Array &); /// get the SolidMechanicsModel::acceleration vector, updated by /// SolidMechanicsModel::updateAcceleration AKANTU_GET_MACRO(Acceleration, *acceleration, Array &); /// get the SolidMechanicsModel::external_force vector (external forces) AKANTU_GET_MACRO(ExternalForce, *external_force, Array &); /// get the SolidMechanicsModel::force vector (external forces) Array & getForce() { AKANTU_DEBUG_WARNING("getForce was maintained for backward compatibility, " "use getExternalForce instead"); return *external_force; } /// get the SolidMechanicsModel::internal_force vector (internal forces) AKANTU_GET_MACRO(InternalForce, *internal_force, Array &); /// get the SolidMechanicsModel::blocked_dofs vector AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); - /// get the value of the SolidMechanicsModel::increment_flag - AKANTU_GET_MACRO(IncrementFlag, increment_flag, bool); - /// get a particular material (by material index) inline Material & getMaterial(UInt mat_index); /// get a particular material (by material index) inline const Material & getMaterial(UInt mat_index) const; /// get a particular material (by material name) inline Material & getMaterial(const std::string & name); /// get a particular material (by material name) inline const Material & getMaterial(const std::string & name) const; /// get a particular material id from is name inline UInt getMaterialIndex(const std::string & name) const; /// give the number of materials inline UInt getNbMaterials() const { return materials.size(); } /// give the material internal index from its id Int getInternalIndexFromID(const ID & id) const; /// compute the stable time step Real getStableTimeStep(); /// get the energies Real getEnergy(const std::string & energy_id); /// compute the energy for energy Real getEnergy(const std::string & energy_id, const ElementType & type, UInt index); AKANTU_GET_MACRO(MaterialByElement, material_index, const ElementTypeMapArray &); AKANTU_GET_MACRO(MaterialLocalNumbering, material_local_numbering, const ElementTypeMapArray &); /// vectors containing local material element index for each global element /// index AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialByElement, material_index, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialByElement, material_index, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialLocalNumbering, material_local_numbering, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialLocalNumbering, material_local_numbering, UInt); AKANTU_GET_MACRO_NOT_CONST(MaterialSelector, *material_selector, MaterialSelector &); AKANTU_SET_MACRO(MaterialSelector, material_selector, std::shared_ptr); /// Access the non_local_manager interface AKANTU_GET_MACRO(NonLocalManager, *non_local_manager, NonLocalManager &); /// get the FEEngine object to integrate or interpolate on the boundary FEEngine & getFEEngineBoundary(const ID & name = "") override; protected: friend class Material; protected: /// compute the stable time step Real getStableTimeStep(const GhostType & ghost_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// conversion coefficient form force/mass to acceleration Real f_m2a; /// displacements array Array * displacement; UInt displacement_release{0}; /// displacements array at the previous time step (used in finite deformation) Array * previous_displacement{nullptr}; /// increment of displacement Array * displacement_increment{nullptr}; /// lumped mass array Array * mass{nullptr}; /// Check if materials need to recompute the mass array bool need_to_reassemble_lumped_mass{true}; /// Check if materials need to recompute the mass matrix bool need_to_reassemble_mass{true}; /// velocities array Array * velocity{nullptr}; /// accelerations array Array * acceleration{nullptr}; /// accelerations array // Array * increment_acceleration; /// external forces array Array * external_force{nullptr}; /// internal forces array Array * internal_force{nullptr}; /// array specifing if a degree of freedom is blocked or not Array * blocked_dofs{nullptr}; /// array of current position used during update residual Array * current_position{nullptr}; UInt current_position_release{0}; /// Arrays containing the material index for each element ElementTypeMapArray material_index; /// Arrays containing the position in the element filter of the material /// (material's local numbering) ElementTypeMapArray material_local_numbering; /// list of used materials std::vector> materials; /// mapping between material name and material internal id std::map materials_names_to_id; /// class defining of to choose a material std::shared_ptr material_selector; /// flag defining if the increment must be computed or not bool increment_flag; /// tells if the material are instantiated bool are_materials_instantiated; using flatten_internal_map = std::map, ElementTypeMapArray *>; /// map a registered internals to be flattened for dump purposes flatten_internal_map registered_internals; /// non local manager std::unique_ptr non_local_manager; }; /* -------------------------------------------------------------------------- */ namespace BC { namespace Neumann { using FromStress = FromHigherDim; using FromTraction = FromSameDim; } // namespace Neumann } // namespace BC } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "parser.hh" #include "solid_mechanics_model_inline_impl.cc" #include "solid_mechanics_model_tmpl.hh" /* -------------------------------------------------------------------------- */ #endif /* __AKANTU_SOLID_MECHANICS_MODEL_HH__ */ diff --git a/src/solver/solver_vector.hh b/src/solver/solver_vector.hh index e783a94d1..53e19a2e4 100644 --- a/src/solver/solver_vector.hh +++ b/src/solver/solver_vector.hh @@ -1,80 +1,87 @@ /** * @file solver_vector.hh * * @author Nicolas Richart * * @date creation Tue Jan 01 2019 * * @brief A Documented file. * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_array.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLVER_VECTOR_HH__ #define __AKANTU_SOLVER_VECTOR_HH__ namespace akantu { class DOFManager; } namespace akantu { class SolverVector { public: SolverVector(DOFManager & dof_manager, const ID & id = "solver_vector") : id(id), _dof_manager(dof_manager) {} SolverVector(const SolverVector & vector, const ID & id = "solver_vector") : id(id), _dof_manager(vector._dof_manager) {} virtual ~SolverVector() = default; // resize the vector to the size of the problem virtual void resize() = 0; // clear the vector virtual void clear() = 0; virtual operator const Array &() const = 0; - virtual Int size() = 0; - virtual Int localSize() = 0; + virtual Int size() const = 0; + virtual Int localSize() const = 0; virtual SolverVector & operator+(const SolverVector & y) = 0; virtual SolverVector & operator=(const SolverVector & y) = 0; UInt & release() { return release_; } UInt release() const { return release_; } + virtual void printself(std::ostream & stream, int indent = 0) const = 0; + protected: ID id; /// Underlying dof manager DOFManager & _dof_manager; UInt release_{0}; }; +inline std::ostream & operator<<(std::ostream & stream, SolverVector & _this) { + _this.printself(stream); + return stream; +} + } // namespace akantu #endif /* __AKANTU_SOLVER_VECTOR_HH__ */ diff --git a/src/solver/solver_vector_default.hh b/src/solver/solver_vector_default.hh index 3367a63f2..7ed5c037a 100644 --- a/src/solver/solver_vector_default.hh +++ b/src/solver/solver_vector_default.hh @@ -1,135 +1,140 @@ /** * @file solver_vector_default.hh * * @author Nicolas Richart * * @date creation Tue Jan 01 2019 * * @brief A Documented file. * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "solver_vector.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLVER_VECTOR_DEFAULT_HH__ #define __AKANTU_SOLVER_VECTOR_DEFAULT_HH__ namespace akantu { class DOFManagerDefault; } // namespace akantu namespace akantu { class SolverVectorArray : public SolverVector { public: SolverVectorArray(DOFManagerDefault & dof_manager, const ID & id); SolverVectorArray(const SolverVectorArray & vector, const ID & id); virtual ~SolverVectorArray() = default; virtual Array & getVector() = 0; virtual const Array & getVector() const = 0; - virtual Int size() = 0; - virtual Int localSize() = 0; + void printself(std::ostream & stream, int indent = 0) const override{ + std::string space(indent, AKANTU_INDENT); + stream << space << "SolverVectorArray [" << std::endl; + stream << space << " + id: " << id << std::endl; + this->getVector().printself(stream, indent + 1); + stream << space << "]" << std::endl; + } }; /* -------------------------------------------------------------------------- */ template class SolverVectorArrayTmpl : public SolverVectorArray { public: SolverVectorArrayTmpl(DOFManagerDefault & dof_manager, Array_ & vector, const ID & id = "solver_vector_default") : SolverVectorArray(dof_manager, id), dof_manager(dof_manager), vector(vector) {} template ::value> * = nullptr> SolverVectorArrayTmpl(DOFManagerDefault & dof_manager, const ID & id = "solver_vector_default") : SolverVectorArray(dof_manager, id), dof_manager(dof_manager), vector(0, 1, id + ":vector") {} SolverVectorArrayTmpl(const SolverVectorArrayTmpl & vector, const ID & id = "solver_vector_default") : SolverVectorArray(vector, id), dof_manager(vector.dof_manager), vector(vector.vector) {} operator const Array &() const override { return getVector(); }; virtual operator Array &() { return getVector(); }; SolverVector & operator+(const SolverVector & y) override; SolverVector & operator=(const SolverVector & y) override; void resize() { static_assert(not std::is_const>::value, "Cannot resize a const Array"); this->vector.resize(this->localSize(), 0.); ++this->release_; } void clear() { static_assert(not std::is_const>::value, "Cannot clear a const Array"); this->vector.clear(); ++this->release_; } public: Array & getVector() override { return vector; } const Array & getVector() const override { return vector; } - Int size() override; - Int localSize() override; + Int size() const override; + Int localSize() const override; virtual Array & getGlobalVector() { return this->vector; } virtual void setGlobalVector(const Array & solution) { this->vector.copy(solution); } protected: DOFManagerDefault & dof_manager; Array_ vector; template friend class SolverVectorArrayTmpl; }; /* -------------------------------------------------------------------------- */ using SolverVectorDefault = SolverVectorArrayTmpl>; /* -------------------------------------------------------------------------- */ template using SolverVectorDefaultWrap = SolverVectorArrayTmpl; template decltype(auto) make_solver_vector_default_wrap(DOFManagerDefault & dof_manager, Array & vector) { return SolverVectorDefaultWrap(dof_manager, vector); } } // namespace akantu /* -------------------------------------------------------------------------- */ #include "solver_vector_default_tmpl.hh" /* -------------------------------------------------------------------------- */ #endif /* __AKANTU_SOLVER_VECTOR_DEFAULT_HH__ */ diff --git a/src/solver/solver_vector_default_tmpl.hh b/src/solver/solver_vector_default_tmpl.hh index 3cff162a0..329d9c4d0 100644 --- a/src/solver/solver_vector_default_tmpl.hh +++ b/src/solver/solver_vector_default_tmpl.hh @@ -1,81 +1,81 @@ /** * @file solver_vector_default_tmpl.hh * * @author Nicolas Richart * * @date creation Tue Jan 01 2019 * * @brief A Documented file. * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "dof_manager_default.hh" #include "solver_vector_default.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLVER_VECTOR_DEFAULT_TMPL_HH__ #define __AKANTU_SOLVER_VECTOR_DEFAULT_TMPL_HH__ namespace akantu { /* -------------------------------------------------------------------------- */ inline SolverVectorArray::SolverVectorArray(DOFManagerDefault & dof_manager, const ID & id) : SolverVector(dof_manager, id) {} /* -------------------------------------------------------------------------- */ inline SolverVectorArray::SolverVectorArray(const SolverVectorArray & vector, const ID & id) : SolverVector(vector, id) {} /* -------------------------------------------------------------------------- */ template SolverVector & SolverVectorArrayTmpl:: operator+(const SolverVector & y) { const auto & y_ = dynamic_cast(y); this->vector += y_.getVector(); ++this->release_; return *this; } /* -------------------------------------------------------------------------- */ template SolverVector & SolverVectorArrayTmpl:: operator=(const SolverVector & y) { const auto & y_ = dynamic_cast(y); this->vector.copy(y_.getVector()); this->release_ = y.release(); return *this; } /* -------------------------------------------------------------------------- */ -template inline Int SolverVectorArrayTmpl::size() { +template inline Int SolverVectorArrayTmpl::size() const { return this->dof_manager.getSystemSize(); } /* -------------------------------------------------------------------------- */ -template inline Int SolverVectorArrayTmpl::localSize() { +template inline Int SolverVectorArrayTmpl::localSize() const { return dof_manager.getLocalSystemSize(); } } // namespace akantu #endif /* __AKANTU_SOLVER_VECTOR_DEFAULT_TMPL_HH__ */ diff --git a/src/solver/solver_vector_petsc.cc b/src/solver/solver_vector_petsc.cc index 0fac5df75..af6ae8a8f 100644 --- a/src/solver/solver_vector_petsc.cc +++ b/src/solver/solver_vector_petsc.cc @@ -1,270 +1,287 @@ /** * @file solver_vector_petsc.cc * * @author Nicolas Richart * * @date creation Tue Jan 01 2019 * * @brief A Documented file. * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "solver_vector_petsc.hh" #include "dof_manager_petsc.hh" #include "mpi_communicator_data.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ SolverVectorPETSc::SolverVectorPETSc(DOFManagerPETSc & dof_manager, const ID & id) : SolverVector(dof_manager, id), dof_manager(dof_manager) { auto && mpi_comm = dof_manager.getMPIComm(); PETSc_call(VecCreate, mpi_comm, &x); detail::PETScSetName(x, id); PETSc_call(VecSetFromOptions, x); auto local_system_size = dof_manager.getLocalSystemSize(); auto nb_local_dofs = dof_manager.getPureLocalSystemSize(); PETSc_call(VecSetSizes, x, nb_local_dofs, PETSC_DECIDE); VecType vec_type; PETSc_call(VecGetType, x, &vec_type); if (std::string(vec_type) == std::string(VECMPI)) { PetscInt lowest_gidx, highest_gidx; PETSc_call(VecGetOwnershipRange, x, &lowest_gidx, &highest_gidx); std::vector ghost_idx; for (auto && d : arange(local_system_size)) { int gidx = dof_manager.localToGlobalEquationNumber(d); if (gidx != -1) { if ((gidx < lowest_gidx) or (gidx >= highest_gidx)) { ghost_idx.push_back(gidx); } } } PETSc_call(VecMPISetGhost, x, ghost_idx.size(), ghost_idx.data()); } else { std::vector idx(nb_local_dofs); std::iota(idx.begin(), idx.end(), 0); ISLocalToGlobalMapping is; PETSc_call(ISLocalToGlobalMappingCreate, PETSC_COMM_SELF, 1, idx.size(), idx.data(), PETSC_COPY_VALUES, &is); PETSc_call(VecSetLocalToGlobalMapping, x, is); PETSc_call(ISLocalToGlobalMappingDestroy, &is); } } /* -------------------------------------------------------------------------- */ SolverVectorPETSc::SolverVectorPETSc(const SolverVectorPETSc & vector, const ID & id) : SolverVector(vector, id), dof_manager(vector.dof_manager) { if (vector.x) { PETSc_call(VecDuplicate, vector.x, &x); PETSc_call(VecCopy, vector.x, x); detail::PETScSetName(x, id); } } +/* -------------------------------------------------------------------------- */ +void SolverVectorPETSc::printself(std::ostream & stream, int indent) const { + std::string space(indent, AKANTU_INDENT); + stream << space << "SolverVectorPETSc [" << std::endl; + stream << space << " + id: " << id << std::endl; + PETSc_call(PetscViewerPushFormat, PETSC_VIEWER_STDOUT_WORLD, + PETSC_VIEWER_ASCII_INDEX); + PETSc_call(VecView, x, PETSC_VIEWER_STDOUT_WORLD); + PETSc_call(PetscViewerPopFormat, PETSC_VIEWER_STDOUT_WORLD); + stream << space << "]" << std::endl; +} + /* -------------------------------------------------------------------------- */ SolverVectorPETSc::SolverVectorPETSc(Vec x, DOFManagerPETSc & dof_manager, const ID & id) : SolverVector(dof_manager, id), dof_manager(dof_manager) { PETSc_call(VecDuplicate, x, &this->x); PETSc_call(VecCopy, x, this->x); detail::PETScSetName(x, id); } /* -------------------------------------------------------------------------- */ SolverVectorPETSc::~SolverVectorPETSc() { if (x) { PETSc_call(VecDestroy, &x); } } /* -------------------------------------------------------------------------- */ void SolverVectorPETSc::resize() { // the arrays are destroyed and recreated in the dof manager // resize is so not implemented } /* -------------------------------------------------------------------------- */ void SolverVectorPETSc::clear() { PETSc_call(VecSet, x, 0.); applyModifications(); } /* -------------------------------------------------------------------------- */ void SolverVectorPETSc::applyModifications() { PETSc_call(VecAssemblyBegin, x); PETSc_call(VecAssemblyEnd, x); updateGhost(); } /* -------------------------------------------------------------------------- */ void SolverVectorPETSc::updateGhost() { Vec x_ghosted{nullptr}; PETSc_call(VecGhostGetLocalForm, x, &x_ghosted); if (x_ghosted) { PETSc_call(VecGhostUpdateBegin, x, INSERT_VALUES, SCATTER_FORWARD); PETSc_call(VecGhostUpdateEnd, x, INSERT_VALUES, SCATTER_FORWARD); } PETSc_call(VecGhostRestoreLocalForm, x, &x_ghosted); } /* -------------------------------------------------------------------------- */ void SolverVectorPETSc::getValues(const Array & idx, Array & values) const { if (idx.size() == 0) return; ISLocalToGlobalMapping is_ltog_map; PETSc_call(VecGetLocalToGlobalMapping, x, &is_ltog_map); PetscInt n; Array lidx(idx.size()); PETSc_call(ISGlobalToLocalMappingApply, is_ltog_map, IS_GTOLM_MASK, idx.size(), idx.storage(), &n, lidx.storage()); getValuesLocal(lidx, values); } /* -------------------------------------------------------------------------- */ void SolverVectorPETSc::getValuesLocal(const Array & idx, Array & values) const { if (idx.size() == 0) return; Vec x_ghosted{nullptr}; PETSc_call(VecGhostGetLocalForm, x, &x_ghosted); // VecScatterBegin(scatter, x, x_local, INSERT_VALUES, SCATTER_FORWARD); // VecScatterEnd(scatter, x, x_local, INSERT_VALUES, SCATTER_FORWARD); if (not x_ghosted) { const PetscScalar * array; PETSc_call(VecGetArrayRead, x, &array); for (auto && data : zip(idx, make_view(values))) { auto i = std::get<0>(data); if (i != -1) { std::get<1>(data) = array[i]; } } PETSc_call(VecRestoreArrayRead, x, &array); return; } PETSc_call(VecSetOption, x_ghosted, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE); PETSc_call(VecGetValues, x_ghosted, idx.size(), idx.storage(), values.storage()); PETSc_call(VecGhostRestoreLocalForm, x, &x_ghosted); } /* -------------------------------------------------------------------------- */ void SolverVectorPETSc::addValues(const Array & gidx, const Array & values, Real scale_factor) { Real * to_add = values.storage(); Array scaled_array; if (scale_factor != 1.) { - scaled_array = values; + scaled_array.copy(values, false); scaled_array *= scale_factor; to_add = scaled_array.storage(); } PETSc_call(VecSetOption, x, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE); PETSc_call(VecSetValues, x, gidx.size(), gidx.storage(), to_add, ADD_VALUES); applyModifications(); } /* -------------------------------------------------------------------------- */ void SolverVectorPETSc::addValuesLocal(const Array & lidx, const Array & values, Real scale_factor) { Vec x_ghosted{nullptr}; PETSc_call(VecGhostGetLocalForm, x, &x_ghosted); if (not x_ghosted) { Real * to_add = values.storage(); Array scaled_array; if (scale_factor != 1.) { - scaled_array = values; + scaled_array.copy(values, false); scaled_array *= scale_factor; to_add = scaled_array.storage(); } - PETSc_call(VecSetOption, x, VEC_IGNORE_NEGATIVE_INDICES, - PETSC_TRUE); + PETSc_call(VecSetOption, x, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE); PETSc_call(VecSetValuesLocal, x, lidx.size(), lidx.storage(), to_add, ADD_VALUES); return; } PETSc_call(VecGhostRestoreLocalForm, x, &x_ghosted); ISLocalToGlobalMapping is_ltog_map; PETSc_call(VecGetLocalToGlobalMapping, x, &is_ltog_map); Array gidx(lidx.size()); PETSc_call(ISLocalToGlobalMappingApply, is_ltog_map, lidx.size(), lidx.storage(), gidx.storage()); addValues(gidx, values, scale_factor); } /* -------------------------------------------------------------------------- */ SolverVectorPETSc::operator const Array &() const { const_cast &>(cache).resize(local_size()); auto xl = internal::make_petsc_local_vector(x); auto cachep = internal::make_petsc_wraped_vector(this->cache); PETSc_call(VecCopy, cachep, xl); return cache; } /* -------------------------------------------------------------------------- */ -SolverVector & SolverVectorPETSc::operator=(const SolverVector & y) { - auto & y_ = dynamic_cast(y); - PETSc_call(VecCopy, y_.x, x); - release_ = y_.release_; +SolverVectorPETSc & SolverVectorPETSc::operator=(const SolverVectorPETSc & y) { + if (size() != y.size()) { + PETSc_call(VecDuplicate, y, &x); + } + + PETSc_call(VecCopy, y.x, x); + release_ = y.release_; return *this; } +/* -------------------------------------------------------------------------- */ +SolverVector & SolverVectorPETSc::operator=(const SolverVector & y) { + const auto & y_ = aka::as_type(y); + return operator=(y_); +} /* -------------------------------------------------------------------------- */ SolverVector & SolverVectorPETSc::operator+(const SolverVector & y) { auto & y_ = dynamic_cast(y); PETSc_call(VecAXPY, x, 1., y_.x); release_ = y_.release_; return *this; } - - } // namespace akantu diff --git a/src/solver/solver_vector_petsc.hh b/src/solver/solver_vector_petsc.hh index 8b32d21b5..b18f302de 100644 --- a/src/solver/solver_vector_petsc.hh +++ b/src/solver/solver_vector_petsc.hh @@ -1,202 +1,205 @@ /** * @file solver_vector_petsc.hh * * @author Nicolas Richart * * @date creation Tue Jan 01 2019 * * @brief A Documented file. * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "dof_manager_petsc.hh" #include "solver_vector.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLVER_VECTOR_PETSC_HH__ #define __AKANTU_SOLVER_VECTOR_PETSC_HH__ namespace akantu { class DOFManagerPETSc; } // namespace akantu namespace akantu { /* -------------------------------------------------------------------------- */ namespace internal { /* ------------------------------------------------------------------------ */ class PETScVector { public: virtual ~PETScVector() = default; operator Vec&() { return x; } operator const Vec&() const { return x; } Int size() const { PetscInt n; PETSc_call(VecGetSize, x, &n); return n; } Int local_size() const { PetscInt n; PETSc_call(VecGetLocalSize, x, &n); return n; } AKANTU_GET_MACRO_NOT_CONST(Vec, x, auto &); AKANTU_GET_MACRO(Vec, x, const auto &); protected: Vec x{nullptr}; }; } // namespace internal /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ class SolverVectorPETSc : public SolverVector, public internal::PETScVector { public: SolverVectorPETSc(DOFManagerPETSc & dof_manager, const ID & id = "solver_vector_petsc"); SolverVectorPETSc(const SolverVectorPETSc & vector, const ID & id = "solver_vector_petsc"); SolverVectorPETSc(Vec vec, DOFManagerPETSc & dof_manager, const ID & id = "solver_vector_petsc"); ~SolverVectorPETSc() override; // resize the vector to the size of the problem void resize() override; void clear() override; operator const Array &() const override; SolverVector & operator+(const SolverVector & y) override; SolverVector & operator=(const SolverVector & y) override; + SolverVectorPETSc & operator=(const SolverVectorPETSc & y); /// get values using processors global indexes void getValues(const Array & idx, Array & values) const; /// get values using processors local indexes void getValuesLocal(const Array & idx, Array & values) const; /// adding values to the vector using the global indices void addValues(const Array & gidx, const Array & values, Real scale_factor = 1.); /// adding values to the vector using the local indices void addValuesLocal(const Array & lidx, const Array & values, Real scale_factor = 1.); - Int size() override { return internal::PETScVector::size(); } - Int localSize() override { return internal::PETScVector::local_size(); } + Int size() const override { return internal::PETScVector::size(); } + Int localSize() const override { return internal::PETScVector::local_size(); } + void printself(std::ostream & stream, int indent = 0) const override; + protected: void applyModifications(); void updateGhost(); protected: DOFManagerPETSc & dof_manager; // used for the conversion operator Array cache; }; /* -------------------------------------------------------------------------- */ namespace internal { /* ------------------------------------------------------------------------ */ template class PETScWrapedVector : public PETScVector { public: PETScWrapedVector(Array && array) : array(array) { PETSc_call(VecCreateSeqWithArray, PETSC_COMM_SELF, 1, array.size(), array.storage(), &x); } ~PETScWrapedVector() { PETSc_call(VecDestroy, &x); } private: Array array; }; /* ------------------------------------------------------------------------ */ template class PETScLocalVector : public PETScVector { public: PETScLocalVector(const Vec & g) : g(g) { PETSc_call(VecGetLocalVectorRead, g, x); } PETScLocalVector(const SolverVectorPETSc & g) : PETScLocalVector(g.getVec()) {} ~PETScLocalVector() { PETSc_call(VecRestoreLocalVectorRead, g, x); PETSc_call(VecDestroy, &x); } private: const Vec & g; }; template <> class PETScLocalVector : public PETScVector { public: PETScLocalVector(Vec & g) : g(g) { PETSc_call(VecGetLocalVectorRead, g, x); } PETScLocalVector(SolverVectorPETSc & g) : PETScLocalVector(g.getVec()) {} ~PETScLocalVector() { PETSc_call(VecRestoreLocalVectorRead, g, x); PETSc_call(VecDestroy, &x); } private: Vec & g; }; /* ------------------------------------------------------------------------ */ template decltype(auto) make_petsc_wraped_vector(Array && array) { return PETScWrapedVector(std::forward(array)); } template < typename V, std::enable_if_t>::value> * = nullptr> decltype(auto) make_petsc_local_vector(V && vec) { constexpr auto read_only = std::is_const>::value; return PETScLocalVector(vec); } template >::value> * = nullptr> decltype(auto) make_petsc_local_vector(V && vec) { constexpr auto read_only = std::is_const>::value; return PETScLocalVector( dynamic_cast &>(vec)); } } // namespace internal } // namespace akantu #endif /* __AKANTU_SOLVER_VECTOR_PETSC_HH__ */ diff --git a/src/solver/sparse_matrix_aij.cc b/src/solver/sparse_matrix_aij.cc index 1711145e8..38a74cd52 100644 --- a/src/solver/sparse_matrix_aij.cc +++ b/src/solver/sparse_matrix_aij.cc @@ -1,290 +1,290 @@ /** * @file sparse_matrix_aij.cc * * @author Nicolas Richart * * @date creation: Fri Aug 21 2015 * @date last modification: Mon Dec 04 2017 * * @brief Implementation of the AIJ sparse matrix * * @section LICENSE * * Copyright (©) 2015-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 "sparse_matrix_aij.hh" #include "aka_iterators.hh" #include "dof_manager_default.hh" #include "dof_synchronizer.hh" #include "solver_vector_default.hh" #include "terms_to_assemble.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ SparseMatrixAIJ::SparseMatrixAIJ(DOFManagerDefault & dof_manager, const MatrixType & matrix_type, const ID & id) : SparseMatrix(dof_manager, matrix_type, id), dof_manager(dof_manager), irn(0, 1, id + ":irn"), jcn(0, 1, id + ":jcn"), a(0, 1, id + ":a") {} /* -------------------------------------------------------------------------- */ SparseMatrixAIJ::SparseMatrixAIJ(const SparseMatrixAIJ & matrix, const ID & id) : SparseMatrix(matrix, id), dof_manager(matrix.dof_manager), irn(matrix.irn, id + ":irn"), jcn(matrix.jcn, id + ":jcn"), a(matrix.a, id + ":a") {} /* -------------------------------------------------------------------------- */ SparseMatrixAIJ::~SparseMatrixAIJ() = default; /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::applyBoundary(Real block_val) { AKANTU_DEBUG_IN(); const auto & blocked_dofs = this->dof_manager.getGlobalBlockedDOFs(); auto begin = blocked_dofs.begin(); auto end = blocked_dofs.end(); auto is_blocked = [&](auto && i) -> bool { auto il = this->dof_manager.globalToLocalEquationNumber(i); return std::binary_search(begin, end, il); }; for (auto && ij_a : zip(irn, jcn, a)) { UInt ni = std::get<0>(ij_a) - 1; UInt nj = std::get<1>(ij_a) - 1; if (is_blocked(ni) or is_blocked(nj)) { // clang-format off std::get<2>(ij_a) = std::get<0>(ij_a) != std::get<1>(ij_a) ? 0. : this->dof_manager.isLocalOrMasterDOF(ni) ? block_val : 0.; // clang-format on } } this->value_release++; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::saveProfile(const std::string & filename) const { AKANTU_DEBUG_IN(); std::ofstream outfile; outfile.open(filename.c_str()); UInt m = this->size_; auto & comm = dof_manager.getCommunicator(); // write header if (comm.whoAmI() == 0) { outfile << "%%MatrixMarket matrix coordinate pattern"; if (this->matrix_type == _symmetric) outfile << " symmetric"; else outfile << " general"; outfile << std::endl; outfile << m << " " << m << " " << this->nb_non_zero << std::endl; } for (auto p : arange(comm.getNbProc())) { // write content if (comm.whoAmI() == p) { for (UInt i = 0; i < this->nb_non_zero; ++i) { outfile << this->irn.storage()[i] << " " << this->jcn.storage()[i] << " 1" << std::endl; } } comm.barrier(); } outfile.close(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::saveMatrix(const std::string & filename) const { AKANTU_DEBUG_IN(); auto & comm = dof_manager.getCommunicator(); // open and set the properties of the stream std::ofstream outfile; if (0 == comm.whoAmI()) { outfile.open(filename.c_str()); } else { outfile.open(filename.c_str(), std::ios_base::app); } outfile.precision(std::numeric_limits::digits10); // write header decltype(nb_non_zero) nnz = this->nb_non_zero; comm.allReduce(nnz); if (comm.whoAmI() == 0) { outfile << "%%MatrixMarket matrix coordinate real"; if (this->matrix_type == _symmetric) outfile << " symmetric"; else outfile << " general"; outfile << std::endl; outfile << this->size_ << " " << this->size_ << " " << nnz << std::endl; } for (auto p : arange(comm.getNbProc())) { // write content if (comm.whoAmI() == p) { for (UInt i = 0; i < this->nb_non_zero; ++i) { outfile << this->irn(i) << " " << this->jcn(i) << " " << this->a(i) << std::endl; } } comm.barrier(); } // time to end outfile.close(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::matVecMul(const Array & x, Array & y, Real alpha, Real beta) const { AKANTU_DEBUG_IN(); y *= beta; auto i_it = this->irn.begin(); auto j_it = this->jcn.begin(); auto a_it = this->a.begin(); auto a_end = this->a.end(); auto x_it = x.begin_reinterpret(x.size() * x.getNbComponent()); auto y_it = y.begin_reinterpret(x.size() * x.getNbComponent()); for (; a_it != a_end; ++i_it, ++j_it, ++a_it) { Int i = this->dof_manager.globalToLocalEquationNumber(*i_it - 1); Int j = this->dof_manager.globalToLocalEquationNumber(*j_it - 1); const Real & A = *a_it; y_it[i] += alpha * A * x_it[j]; if ((this->matrix_type == _symmetric) && (i != j)) y_it[j] += alpha * A * x_it[i]; } if (this->dof_manager.hasSynchronizer()) this->dof_manager.getSynchronizer().reduceSynchronize(y); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::matVecMul(const SolverVector & _x, SolverVector & _y, Real alpha, Real beta) const { AKANTU_DEBUG_IN(); - auto && x = dynamic_cast(_x).getVector(); - auto && y = dynamic_cast(_y).getVector(); + auto && x = aka::as_type(_x).getVector(); + auto && y = aka::as_type(_y).getVector(); this->matVecMul(x, y, alpha, beta); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::copyContent(const SparseMatrix & matrix) { AKANTU_DEBUG_IN(); - const auto & mat = dynamic_cast(matrix); + const auto & mat = aka::as_type(matrix); AKANTU_DEBUG_ASSERT(nb_non_zero == mat.getNbNonZero(), "The to matrix don't have the same profiles"); memcpy(a.storage(), mat.getA().storage(), nb_non_zero * sizeof(Real)); this->value_release++; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::copyProfile(const SparseMatrix & other) { - auto & A = dynamic_cast(other); + auto & A = aka::as_type(other); SparseMatrix::clearProfile(); this->irn = A.irn; this->jcn = A.jcn; this->irn_jcn_k.clear(); UInt i, j, k; for(auto && data : enumerate(irn, jcn)) { std::tie(k, i, j) = data; this->irn_jcn_k[this->key(i - 1, j - 1)] = k; } this->nb_non_zero = this->irn.size(); this->a.resize(this->nb_non_zero); this->a.set(0.); this->size_ = A.size_; this->profile_release++; this->value_release++; } /* -------------------------------------------------------------------------- */ template void SparseMatrixAIJ::addMeToTemplated(MatrixType & B, Real alpha) const { UInt i, j; Real A_ij; for (auto && tuple : zip(irn, jcn, a)) { std::tie(i, j, A_ij) = tuple; B.add(i - 1, j - 1, alpha * A_ij); } } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::addMeTo(SparseMatrix & B, Real alpha) const { if (auto * B_aij = dynamic_cast(&B)) { this->addMeToTemplated(*B_aij, alpha); } else { // this->addMeToTemplated(*this, alpha); } } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::mul(Real alpha) { this->a *= alpha; this->value_release++; } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::clear() { a.set(0.); this->value_release++; } } // namespace akantu diff --git a/src/solver/sparse_matrix_petsc.cc b/src/solver/sparse_matrix_petsc.cc index 5ef0d05d9..513e072ee 100644 --- a/src/solver/sparse_matrix_petsc.cc +++ b/src/solver/sparse_matrix_petsc.cc @@ -1,285 +1,291 @@ /** * @file sparse_matrix_petsc.cc * * @author Aurelia Isabel Cuba Ramos * * @date creation: Mon Dec 13 2010 * @date last modification: Sat Feb 03 2018 * * @brief Implementation of PETSc matrix 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 "sparse_matrix_petsc.hh" #include "dof_manager_petsc.hh" #include "mpi_communicator_data.hh" #include "solver_vector_petsc.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ SparseMatrixPETSc::SparseMatrixPETSc(DOFManagerPETSc & dof_manager, const MatrixType & matrix_type, const ID & id) : SparseMatrix(dof_manager, matrix_type, id), dof_manager(dof_manager) { AKANTU_DEBUG_IN(); auto mpi_comm = dof_manager.getMPIComm(); PETSc_call(MatCreate, mpi_comm, &mat); detail::PETScSetName(mat, id); resize(); PETSc_call(MatSetFromOptions, mat); PETSc_call(MatSetUp, mat); PETSc_call(MatSetOption, mat, MAT_ROW_ORIENTED, PETSC_TRUE); PETSc_call(MatSetOption, mat, MAT_NEW_NONZERO_LOCATIONS, PETSC_TRUE); if (matrix_type == _symmetric) PETSc_call(MatSetOption, mat, MAT_SYMMETRIC, PETSC_TRUE); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SparseMatrixPETSc::SparseMatrixPETSc(const SparseMatrixPETSc & matrix, const ID & id) : SparseMatrix(matrix, id), dof_manager(matrix.dof_manager) { PETSc_call(MatDuplicate, matrix.mat, MAT_COPY_VALUES, &mat); detail::PETScSetName(mat, id); } /* -------------------------------------------------------------------------- */ SparseMatrixPETSc::~SparseMatrixPETSc() { AKANTU_DEBUG_IN(); if (mat) PETSc_call(MatDestroy, &mat); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixPETSc::resize() { auto local_size = dof_manager.getPureLocalSystemSize(); PETSc_call(MatSetSizes, mat, local_size, local_size, size_, size_); auto & is_ltog_mapping = dof_manager.getISLocalToGlobalMapping(); PETSc_call(MatSetLocalToGlobalMapping, mat, is_ltog_mapping, is_ltog_mapping); } /* -------------------------------------------------------------------------- */ /** * Method to save the nonzero pattern and the values stored at each position * @param filename name of the file in which the information will be stored */ void SparseMatrixPETSc::saveMatrix(const std::string & filename) const { AKANTU_DEBUG_IN(); auto mpi_comm = dof_manager.getMPIComm(); /// create Petsc viewer PetscViewer viewer; PETSc_call(PetscViewerASCIIOpen, mpi_comm, filename.c_str(), &viewer); PETSc_call(PetscViewerPushFormat, viewer, PETSC_VIEWER_ASCII_MATRIXMARKET); PETSc_call(MatView, mat, viewer); PETSc_call(PetscViewerPopFormat, viewer); PETSc_call(PetscViewerDestroy, &viewer); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /// Equivalent of *gemv in blas void SparseMatrixPETSc::matVecMul(const SolverVector & _x, SolverVector & _y, Real alpha, Real beta) const { auto & x = aka::as_type(_x); auto & y = aka::as_type(_y); // y = alpha A x + beta y SolverVectorPETSc w(x, this->id + ":tmp"); // w = A x - PETSc_call(MatMult, mat, x.getVec(), w.getVec()); + if (release == 0) { + PETSc_call(VecZeroEntries, w); + } else { + PETSc_call(MatMult, mat, x, w); + } + if (alpha != 1.) { // w = alpha w - PETSc_call(VecScale, w.getVec(), alpha); + PETSc_call(VecScale, w, alpha); } // y = w + beta y - PETSc_call(VecAYPX, y.getVec(), beta, w.getVec()); + PETSc_call(VecAYPX, y, beta, w); } /* -------------------------------------------------------------------------- */ void SparseMatrixPETSc::addMeToImpl(SparseMatrixPETSc & B, Real alpha) const { PETSc_call(MatAXPY, B.mat, alpha, mat, SAME_NONZERO_PATTERN); B.release++; } /* -------------------------------------------------------------------------- */ /** * Method to add another PETSc matrix to this PETSc matrix * @param matrix PETSc matrix to be added * @param alpha the factor specifying how many times the matrix should be added */ void SparseMatrixPETSc::addMeTo(SparseMatrix & B, Real alpha) const { - if (auto * B_petsc = dynamic_cast(&B)) { - this->addMeToImpl(*B_petsc, alpha); + if (aka::is_of_type(B)) { + auto & B_petsc = aka::as_type(B); + this->addMeToImpl(B_petsc, alpha); } else { AKANTU_TO_IMPLEMENT(); // this->addMeToTemplated(*this, alpha); } } /* -------------------------------------------------------------------------- */ /** * MatSetValues() generally caches the values. The matrix is ready to * use only after MatAssemblyBegin() and MatAssemblyEnd() have been * called. (http://www.mcs.anl.gov/petsc/) */ void SparseMatrixPETSc::applyModifications() { this->beginAssembly(); this->endAssembly(); } /* -------------------------------------------------------------------------- */ void SparseMatrixPETSc::beginAssembly() { PETSc_call(MatAssemblyBegin, mat, MAT_FINAL_ASSEMBLY); } /* -------------------------------------------------------------------------- */ void SparseMatrixPETSc::endAssembly() { PETSc_call(MatAssemblyEnd, mat, MAT_FINAL_ASSEMBLY); PETSc_call(MatSetOption, mat, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE); this->release++; } /* -------------------------------------------------------------------------- */ void SparseMatrixPETSc::copyProfile(const SparseMatrix & other) { - auto & A = dynamic_cast(other); + auto & A = aka::as_type(other); MatDestroy(&mat); MatDuplicate(A.mat, MAT_DO_NOT_COPY_VALUES, &mat); } /* -------------------------------------------------------------------------- */ void SparseMatrixPETSc::applyBoundary(Real block_val) { AKANTU_DEBUG_IN(); const auto & blocked_dofs = this->dof_manager.getGlobalBlockedDOFs(); // std::vector rows; // for (auto && data : enumerate(blocked)) { // if (std::get<1>(data)) { // rows.push_back(std::get<0>(data)); // } // } //applyModifications(); static int c = 0; saveMatrix("before_blocked_" + std::to_string(c) + ".mtx"); PETSc_call(MatZeroRowsColumnsLocal, mat, blocked_dofs.size(), blocked_dofs.storage(), block_val, nullptr, nullptr); saveMatrix("after_blocked_" + std::to_string(c) + ".mtx"); ++c; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixPETSc::mul(Real alpha) { PETSc_call(MatScale, mat, alpha); this->release++; } /* -------------------------------------------------------------------------- */ void SparseMatrixPETSc::clear() { PETSc_call(MatZeroEntries, mat); this->release++; } /* -------------------------------------------------------------------------- */ void SparseMatrixPETSc::clearProfile() { SparseMatrix::clearProfile(); PETSc_call(MatResetPreallocation, mat); PETSc_call(MatSetOption, mat, MAT_NEW_NONZERO_LOCATIONS, PETSC_TRUE); // PETSc_call(MatSetOption, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE); // PETSc_call(MatSetOption, MAT_NEW_NONZERO_ALLOCATIONS, PETSC_TRUE); // PETSc_call(MatSetOption, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE); clear(); } /* -------------------------------------------------------------------------- */ UInt SparseMatrixPETSc::add(UInt i, UInt j) { PETSc_call(MatSetValue, mat, i, j, 0, ADD_VALUES); return 0; } /* -------------------------------------------------------------------------- */ void SparseMatrixPETSc::add(UInt i, UInt j, Real val) { PETSc_call(MatSetValue, mat, i, j, val, ADD_VALUES); } /* -------------------------------------------------------------------------- */ void SparseMatrixPETSc::addLocal(UInt i, UInt j) { PETSc_call(MatSetValueLocal, mat, i, j, 0, ADD_VALUES); } /* -------------------------------------------------------------------------- */ void SparseMatrixPETSc::addLocal(UInt i, UInt j, Real val) { PETSc_call(MatSetValueLocal, mat, i, j, val, ADD_VALUES); } /* -------------------------------------------------------------------------- */ void SparseMatrixPETSc::addLocal(const Vector & rows, const Vector & cols, const Matrix & vals) { PETSc_call(MatSetValuesLocal, mat, rows.size(), rows.storage(), cols.size(), cols.storage(), vals.storage(), ADD_VALUES); } /* -------------------------------------------------------------------------- */ void SparseMatrixPETSc::addValues(const Vector & rows, const Vector & cols, const Matrix & vals, MatrixType type) { if (type == _unsymmetric and matrix_type == _symmetric) { PETSc_call(MatSetOption, mat, MAT_SYMMETRIC, PETSC_FALSE); PETSc_call(MatSetOption, mat, MAT_STRUCTURALLY_SYMMETRIC, PETSC_FALSE); } PETSc_call(MatSetValues, mat, rows.size(), rows.storage(), cols.size(), cols.storage(), vals.storage(), ADD_VALUES); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/test/test_model/test_common/test_model_solver/test_model_solver_dynamic.cc b/test/test_model/test_common/test_model_solver/test_model_solver_dynamic.cc index 33d6bf58a..feb9a7319 100644 --- a/test/test_model/test_common/test_model_solver/test_model_solver_dynamic.cc +++ b/test/test_model/test_common/test_model_solver/test_model_solver_dynamic.cc @@ -1,239 +1,240 @@ /** * @file test_model_solver_dynamic.cc * * @author Nicolas Richart * * @date creation: Wed Apr 13 2016 * @date last modification: Tue Feb 20 2018 * * @brief Test default dof manager * * @section LICENSE * * Copyright (©) 2016-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 "communicator.hh" #include "element_group.hh" #include "mesh.hh" #include "mesh_accessor.hh" #include "non_linear_solver.hh" /* -------------------------------------------------------------------------- */ #include "dumpable_inline_impl.hh" #include "dumper_element_partition.hh" #include "dumper_iohelper_paraview.hh" /* -------------------------------------------------------------------------- */ #include "test_model_solver_my_model.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef EXPLICIT #define EXPLICIT true #endif using namespace akantu; class Sinusoidal : public BC::Dirichlet::DirichletFunctor { public: Sinusoidal(MyModel & model, Real amplitude, Real pulse_width, Real t) : model(model), A(amplitude), k(2 * M_PI / pulse_width), t(t), v{std::sqrt(model.E / model.rho)} {} void operator()(UInt n, Vector & /*flags*/, Vector & disp, const Vector & coord) const { auto x = coord(_x); model.velocity(n, _x) = k * v * A * sin(k * (x - v * t)); disp(_x) = A * cos(k * (x - v * t)); } private: MyModel & model; Real A{1.}; Real k{2 * M_PI}; Real t{1.}; Real v{1.}; }; static void genMesh(Mesh & mesh, UInt nb_nodes); /* -------------------------------------------------------------------------- */ int main(int argc, char * argv[]) { initialize(argc, argv); UInt prank = Communicator::getStaticCommunicator().whoAmI(); UInt global_nb_nodes = 201; UInt max_steps = 400; Real time_step = 0.001; Mesh mesh(1); Real F = -9.81; bool _explicit = EXPLICIT; const Real pulse_width = 0.2; const Real A = 0.01; if (prank == 0) genMesh(mesh, global_nb_nodes); mesh.distribute(); // mesh.makePeriodic(_x); MyModel model(F, mesh, _explicit); model.forces.clear(); model.blocked.clear(); model.applyBC(Sinusoidal(model, A, pulse_width, 0.), "all"); model.applyBC(BC::Dirichlet::FlagOnly(_x), "border"); if (!_explicit) { model.getNewSolver("dynamic", TimeStepSolverType::_dynamic, NonLinearSolverType::_newton_raphson); model.setIntegrationScheme("dynamic", "disp", IntegrationSchemeType::_trapezoidal_rule_2, IntegrationScheme::_displacement); } else { model.getNewSolver("dynamic", TimeStepSolverType::_dynamic_lumped, NonLinearSolverType::_lumped); model.setIntegrationScheme("dynamic", "disp", IntegrationSchemeType::_central_difference, IntegrationScheme::_acceleration); } model.setTimeStep(time_step); if (prank == 0) { std::cout << std::scientific; std::cout << std::setw(14) << "time" << "," << std::setw(14) << "wext" << "," << std::setw(14) << "epot" << "," << std::setw(14) << "ekin" << "," << std::setw(14) << "total" << "," << std::setw(14) << "max_disp" << "," << std::setw(14) << "min_disp" << std::endl; } Real wext = 0.; model.getDOFManager().clearResidual(); model.assembleResidual(); Real epot = 0; // model.getPotentialEnergy(); Real ekin = 0; // model.getKineticEnergy(); Real einit = ekin + epot; Real etot = ekin + epot - wext - einit; Real max_disp = 0., min_disp = 0.; for (auto && disp : model.displacement) { max_disp = std::max(max_disp, disp); min_disp = std::min(min_disp, disp); } if (prank == 0) { std::cout << std::setw(14) << 0. << "," << std::setw(14) << wext << "," << std::setw(14) << epot << "," << std::setw(14) << ekin << "," << std::setw(14) << etot << "," << std::setw(14) << max_disp << "," << std::setw(14) << min_disp << std::endl; } #if EXPLICIT == false NonLinearSolver & solver = model.getDOFManager().getNonLinearSolver("dynamic"); solver.set("max_iterations", 20); #endif auto * dumper = new DumperParaview("dynamic", "./paraview"); mesh.registerExternalDumper(*dumper, "dynamic", true); mesh.addDumpMesh(mesh); mesh.addDumpFieldExternalToDumper("dynamic", "displacement", model.displacement); mesh.addDumpFieldExternalToDumper("dynamic", "velocity", model.velocity); mesh.addDumpFieldExternalToDumper("dynamic", "forces", model.forces); + mesh.addDumpFieldExternalToDumper("dynamic", "internal_forces", model.internal_forces); mesh.addDumpFieldExternalToDumper("dynamic", "acceleration", model.acceleration); mesh.dump(); for (UInt i = 1; i < max_steps + 1; ++i) { model.applyBC(Sinusoidal(model, A, pulse_width, time_step * (i - 1)), "border"); model.solveStep("dynamic"); mesh.dump(); epot = model.getPotentialEnergy(); ekin = model.getKineticEnergy(); wext += model.getExternalWorkIncrement(); etot = ekin + epot - wext - einit; Real max_disp = 0., min_disp = 0.; for (auto && disp : model.displacement) { max_disp = std::max(max_disp, disp); min_disp = std::min(min_disp, disp); } if (prank == 0) { std::cout << std::setw(14) << time_step * i << "," << std::setw(14) << wext << "," << std::setw(14) << epot << "," << std::setw(14) << ekin << "," << std::setw(14) << etot << "," << std::setw(14) << max_disp << "," << std::setw(14) << min_disp << std::endl; } } // output.close(); finalize(); return EXIT_SUCCESS; } /* -------------------------------------------------------------------------- */ void genMesh(Mesh & mesh, UInt nb_nodes) { MeshAccessor mesh_accessor(mesh); Array & nodes = mesh_accessor.getNodes(); Array & conn = mesh_accessor.getConnectivity(_segment_2); nodes.resize(nb_nodes); auto & all = mesh.createNodeGroup("all_nodes"); for (UInt n = 0; n < nb_nodes; ++n) { nodes(n, _x) = n * (1. / (nb_nodes - 1)); all.add(n); } mesh.createElementGroupFromNodeGroup("all", "all_nodes"); conn.resize(nb_nodes - 1); for (UInt n = 0; n < nb_nodes - 1; ++n) { conn(n, 0) = n; conn(n, 1) = n + 1; } Array & conn_points = mesh_accessor.getConnectivity(_point_1); conn_points.resize(2); conn_points(0, 0) = 0; conn_points(1, 0) = nb_nodes - 1; auto & border = mesh.createElementGroup("border", 0); border.add({_point_1, 0, _not_ghost}, true); border.add({_point_1, 1, _not_ghost}, true); mesh_accessor.makeReady(); } diff --git a/test/test_model/test_common/test_model_solver/test_model_solver_my_model.hh b/test/test_model/test_common/test_model_solver/test_model_solver_my_model.hh index fc33ec22a..0f88c3f59 100644 --- a/test/test_model/test_common/test_model_solver/test_model_solver_my_model.hh +++ b/test/test_model/test_common/test_model_solver/test_model_solver_my_model.hh @@ -1,451 +1,451 @@ /** * @file test_model_solver_my_model.hh * * @author Nicolas Richart * * @date creation: Wed Apr 13 2016 * @date last modification: Tue Feb 20 2018 * * @brief Test default dof manager * * @section LICENSE * * Copyright (©) 2016-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 "aka_iterators.hh" #include "boundary_condition.hh" #include "communicator.hh" #include "data_accessor.hh" #include "dof_manager_default.hh" #include "element_synchronizer.hh" #include "mesh.hh" #include "model_solver.hh" #include "periodic_node_synchronizer.hh" #include "solver_vector_default.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ namespace akantu { #ifndef __AKANTU_TEST_MODEL_SOLVER_MY_MODEL_HH__ #define __AKANTU_TEST_MODEL_SOLVER_MY_MODEL_HH__ /** * =\o-----o-----o-> F * | | * |---- L ----| */ class MyModel : public ModelSolver, public BoundaryCondition, public DataAccessor { public: MyModel(Real F, Mesh & mesh, bool lumped) : ModelSolver(mesh, ModelType::_model, "model_solver", 0), nb_dofs(mesh.getNbNodes()), nb_elements(mesh.getNbElement(_segment_2)), lumped(lumped), E(1.), A(1.), rho(1.), mesh(mesh), displacement(nb_dofs, 1, "disp"), velocity(nb_dofs, 1, "velo"), acceleration(nb_dofs, 1, "accel"), blocked(nb_dofs, 1, "blocked"), forces(nb_dofs, 1, "force_ext"), internal_forces(nb_dofs, 1, "force_int"), stresses(nb_elements, 1, "stress"), strains(nb_elements, 1, "strain"), initial_lengths(nb_elements, 1, "L0") { this->initBC(*this, displacement, forces); - //this->initDOFManager("petsc"); - this->initDOFManager("default"); + this->initDOFManager("petsc"); + //this->initDOFManager("default"); this->getDOFManager().registerDOFs("disp", displacement, _dst_nodal); this->getDOFManager().registerDOFsDerivative("disp", 1, velocity); this->getDOFManager().registerDOFsDerivative("disp", 2, acceleration); this->getDOFManager().registerBlockedDOFs("disp", blocked); displacement.set(0.); velocity.set(0.); acceleration.set(0.); forces.set(0.); blocked.set(false); UInt global_nb_nodes = mesh.getNbGlobalNodes(); for (auto && n : arange(nb_dofs)) { auto global_id = mesh.getNodeGlobalId(n); if (global_id == (global_nb_nodes - 1)) forces(n, _x) = F; if (global_id == 0) blocked(n, _x) = true; } auto cit = this->mesh.getConnectivity(_segment_2).begin(2); auto cend = this->mesh.getConnectivity(_segment_2).end(2); auto L_it = this->initial_lengths.begin(); for (; cit != cend; ++cit, ++L_it) { const Vector & conn = *cit; UInt n1 = conn(0); UInt n2 = conn(1); Real p1 = this->mesh.getNodes()(n1, _x); Real p2 = this->mesh.getNodes()(n2, _x); *L_it = std::abs(p2 - p1); } this->registerDataAccessor(*this); this->registerSynchronizer( const_cast(this->mesh.getElementSynchronizer()), SynchronizationTag::_user_1); } void assembleLumpedMass() { auto & M = this->getDOFManager().getLumpedMatrix("M"); M.clear(); this->assembleLumpedMass(_not_ghost); if (this->mesh.getNbElement(_segment_2, _ghost) > 0) this->assembleLumpedMass(_ghost); is_lumped_mass_assembled = true; } void assembleLumpedMass(const GhostType & ghost_type) { Array M(nb_dofs, 1, 0.); Array m_all_el(this->mesh.getNbElement(_segment_2, ghost_type), 2); for (auto && data : zip(make_view(this->mesh.getConnectivity(_segment_2), 2), make_view(m_all_el, 2))) { const auto & conn = std::get<0>(data); auto & m_el = std::get<1>(data); UInt n1 = conn(0); UInt n2 = conn(1); Real p1 = this->mesh.getNodes()(n1, _x); Real p2 = this->mesh.getNodes()(n2, _x); Real L = std::abs(p2 - p1); Real M_n = rho * A * L / 2; m_el(0) = m_el(1) = M_n; } this->getDOFManager().assembleElementalArrayLocalArray( m_all_el, M, _segment_2, ghost_type); this->getDOFManager().assembleToLumpedMatrix("disp", M, "M"); } void assembleMass() { SparseMatrix & M = this->getDOFManager().getMatrix("M"); M.clear(); Array m_all_el(this->nb_elements, 4); Matrix m(2, 2); m(0, 0) = m(1, 1) = 2; m(0, 1) = m(1, 0) = 1; // under integrated // m(0, 0) = m(1, 1) = 3./2.; // m(0, 1) = m(1, 0) = 3./2.; // lumping the mass matrix // m(0, 0) += m(0, 1); // m(1, 1) += m(1, 0); // m(0, 1) = m(1, 0) = 0; for (auto && data : zip(make_view(this->mesh.getConnectivity(_segment_2), 2), make_view(m_all_el, 2, 2))) { const auto & conn = std::get<0>(data); auto & m_el = std::get<1>(data); UInt n1 = conn(0); UInt n2 = conn(1); Real p1 = this->mesh.getNodes()(n1, _x); Real p2 = this->mesh.getNodes()(n2, _x); Real L = std::abs(p2 - p1); m_el = m; m_el *= rho * A * L / 6.; } this->getDOFManager().assembleElementalMatricesToMatrix( "M", "disp", m_all_el, _segment_2); is_mass_assembled = true; } MatrixType getMatrixType(const ID &) { return _symmetric; } void assembleMatrix(const ID & matrix_id) { if (matrix_id == "K") { if (not is_stiffness_assembled) this->assembleStiffness(); } else if (matrix_id == "M") { if (not is_mass_assembled) this->assembleMass(); } else if (matrix_id == "C") { // pass, no damping matrix } else { AKANTU_EXCEPTION("This solver does not know what to do with a matrix " << matrix_id); } } void assembleLumpedMatrix(const ID & matrix_id) { if (matrix_id == "M") { if (not is_lumped_mass_assembled) this->assembleLumpedMass(); } else { AKANTU_EXCEPTION("This solver does not know what to do with a matrix " << matrix_id); } } void assembleStiffness() { SparseMatrix & K = this->getDOFManager().getMatrix("K"); K.clear(); Matrix k(2, 2); k(0, 0) = k(1, 1) = 1; k(0, 1) = k(1, 0) = -1; Array k_all_el(this->nb_elements, 4); auto k_it = k_all_el.begin(2, 2); auto cit = this->mesh.getConnectivity(_segment_2).begin(2); auto cend = this->mesh.getConnectivity(_segment_2).end(2); for (; cit != cend; ++cit, ++k_it) { const auto & conn = *cit; UInt n1 = conn(0); UInt n2 = conn(1); Real p1 = this->mesh.getNodes()(n1, _x); Real p2 = this->mesh.getNodes()(n2, _x); Real L = std::abs(p2 - p1); auto & k_el = *k_it; k_el = k; k_el *= E * A / L; } this->getDOFManager().assembleElementalMatricesToMatrix( "K", "disp", k_all_el, _segment_2); is_stiffness_assembled = true; } void assembleResidual() { this->getDOFManager().assembleToResidual("disp", forces); internal_forces.clear(); this->assembleResidual(_not_ghost); this->synchronize(SynchronizationTag::_user_1); this->getDOFManager().assembleToResidual("disp", internal_forces, -1.); } void assembleResidual(const GhostType & ghost_type) { Array forces_internal_el( this->mesh.getNbElement(_segment_2, ghost_type), 2); auto cit = this->mesh.getConnectivity(_segment_2, ghost_type).begin(2); auto cend = this->mesh.getConnectivity(_segment_2, ghost_type).end(2); auto f_it = forces_internal_el.begin(2); auto strain_it = this->strains.begin(); auto stress_it = this->stresses.begin(); auto L_it = this->initial_lengths.begin(); for (; cit != cend; ++cit, ++f_it, ++strain_it, ++stress_it, ++L_it) { const auto & conn = *cit; UInt n1 = conn(0); UInt n2 = conn(1); Real u1 = this->displacement(n1, _x); Real u2 = this->displacement(n2, _x); *strain_it = (u2 - u1) / *L_it; *stress_it = E * *strain_it; Real f_n = A * *stress_it; Vector & f = *f_it; f(0) = -f_n; f(1) = f_n; } this->getDOFManager().assembleElementalArrayLocalArray( forces_internal_el, internal_forces, _segment_2, ghost_type); } Real getPotentialEnergy() { Real res = 0; if (not lumped) { res = this->mulVectMatVect(this->displacement, "K", this->displacement); } else { auto strain_it = this->strains.begin(); auto stress_it = this->stresses.begin(); auto strain_end = this->strains.end(); auto L_it = this->initial_lengths.begin(); for (; strain_it != strain_end; ++strain_it, ++stress_it, ++L_it) { res += *strain_it * *stress_it * A * *L_it; } mesh.getCommunicator().allReduce(res, SynchronizerOperation::_sum); } return res / 2.; } Real getKineticEnergy() { Real res = 0; if (not lumped) { res = this->mulVectMatVect(this->velocity, "M", this->velocity); } else { Array & m = dynamic_cast( this->getDOFManager().getLumpedMatrix("M")); auto it = velocity.begin(); auto end = velocity.end(); auto m_it = m.begin(); for (UInt node = 0; it != end; ++it, ++m_it, ++node) { if (mesh.isLocalOrMasterNode(node)) res += *m_it * *it * *it; } mesh.getCommunicator().allReduce(res, SynchronizerOperation::_sum); } return res / 2.; } Real getExternalWorkIncrement() { Real res = 0; auto it = velocity.begin(); auto end = velocity.end(); auto if_it = internal_forces.begin(); auto ef_it = forces.begin(); auto b_it = blocked.begin(); for (UInt node = 0; it != end; ++it, ++if_it, ++ef_it, ++b_it, ++node) { if (mesh.isLocalOrMasterNode(node)) res += (*b_it ? -*if_it : *ef_it) * *it; } mesh.getCommunicator().allReduce(res, SynchronizerOperation::_sum); return res * this->getTimeStep(); } Real mulVectMatVect(const Array & x, const ID & A_id, const Array & y) { Array Ay(nb_dofs); this->getDOFManager().assembleMatMulVectToArray("disp", A_id, y, Ay); Real res = 0.; for (auto && data : zip(arange(nb_dofs), make_view(Ay), make_view(x))) { res += std::get<2>(data) * std::get<1>(data) * mesh.isLocalOrMasterNode(std::get<0>(data)); } mesh.getCommunicator().allReduce(res, SynchronizerOperation::_sum); return res; } void predictor() {} void corrector() {} /* ------------------------------------------------------------------------ */ UInt getNbData(const Array & elements, const SynchronizationTag &) const { return elements.size() * sizeof(Real); } void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const { if (tag == SynchronizationTag::_user_1) { for (const auto & el : elements) { buffer << this->stresses(el.element); } } } void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) { if (tag == SynchronizationTag::_user_1) { auto cit = this->mesh.getConnectivity(_segment_2, _ghost).begin(2); for (const auto & el : elements) { Real stress; buffer >> stress; Real f = A * stress; Vector conn = cit[el.element]; this->internal_forces(conn(0), _x) += -f; this->internal_forces(conn(1), _x) += f; } } } const Mesh & getMesh() const { return mesh; } UInt getSpatialDimension() const { return 1; } auto & getBlockedDOFs() { return blocked; } private: UInt nb_dofs; UInt nb_elements; bool lumped; bool is_stiffness_assembled{false}; bool is_mass_assembled{false}; bool is_lumped_mass_assembled{false}; public: Real E, A, rho; Mesh & mesh; Array displacement; Array velocity; Array acceleration; Array blocked; Array forces; Array internal_forces; Array stresses; Array strains; Array initial_lengths; }; #endif /* __AKANTU_TEST_MODEL_SOLVER_MY_MODEL_HH__ */ } // namespace akantu