diff --git a/cmake/AkantuMacros.cmake b/cmake/AkantuMacros.cmake index 1d0ca8b57..3e784c762 100644 --- a/cmake/AkantuMacros.cmake +++ b/cmake/AkantuMacros.cmake @@ -1,268 +1,273 @@ #=============================================================================== # @file AkantuMacros.cmake # # @author Guillaume Anciaux # @author Nicolas Richart # # @date creation: Fri Oct 22 2010 # @date last modification: Tue Jan 19 2016 # # @brief Set of macros used by akantu cmake files # # @section LICENSE # # Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de # Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des # Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # #=============================================================================== #=============================================================================== function(set_third_party_shared_libirary_name _var _lib) set(${_var} ${PROJECT_BINARY_DIR}/third-party/lib/${CMAKE_SHARED_LIBRARY_PREFIX}${_lib}${CMAKE_SHARED_LIBRARY_SUFFIX} CACHE FILEPATH "" FORCE) endfunction() # ============================================================================== function(_add_file_to_copy target file) get_filename_component(_file_name_we ${file} NAME_WE) get_filename_component(_file_name_ext ${file} EXT) get_filename_component(_file_name ${file} NAME) get_filename_component(_file_path ${file} ABSOLUTE BASE_DIR ${CMAKE_CURRENT_SOURCE_DIR}) - set(copy_target copy_${_file_name_we}_${_file_name_ext}_${target}) - add_custom_target(${copy_target} - DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/${_file_name}) - add_custom_command( - OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/${_file_name} - COMMAND ${CMAKE_COMMAND} -E copy_if_different - ${file} - ${CMAKE_CURRENT_BINARY_DIR} - WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} - DEPENDS ${_file_path} - COMMENT "Copying file ${_file_name} for the target ${target}" - ) - add_dependencies(${target} ${copy_target}) + configure_file( + ${_file_path} + ${CMAKE_CURRENT_BINARY_DIR}/${_file_name} + COPYONLY) + + # set(copy_target copy_${_file_name_we}_${_file_name_ext}_${target}) + # add_custom_target(${copy_target} + # DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/${_file_name}) + # add_custom_command( + # OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/${_file_name} + # COMMAND ${CMAKE_COMMAND} -E copy_if_different + # ${file} + # ${CMAKE_CURRENT_BINARY_DIR} + # WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} + # DEPENDS ${_file_path} + # COMMENT "Copying file ${_file_name} for the target ${target}" + # ) + # add_dependencies(${target} ${copy_target}) endfunction() # ============================================================================== function(get_target_list_of_associated_files tgt files) if(TARGET ${tgt}) get_target_property(_type ${tgt} TYPE) else() set(_type ${tgt}-NOTFOUND) endif() if(_type STREQUAL "SHARED_LIBRARY" OR _type STREQUAL "STATIC_LIBRARY" OR _type STREQUAL "MODULE_LIBRARY" OR _type STREQUAL "EXECUTABLE") get_target_property(_srcs ${tgt} SOURCES) set(_dep_ressources) foreach(_file ${_srcs}) list(APPEND _dep_ressources ${CMAKE_CURRENT_SOURCE_DIR}/${_file}) endforeach() elseif(_type) get_target_property(_dep_ressources ${tgt} RESSOURCES) endif() set(${files} ${_dep_ressources} PARENT_SCOPE) endfunction() #=============================================================================== # Generate the list of currently loaded materials function(generate_material_list) message(STATUS "Determining the list of recognized materials...") package_get_all_include_directories( AKANTU_LIBRARY_INCLUDE_DIRS ) package_get_all_external_informations( PRIVATE_INCLUDE AKANTU_PRIVATE_EXTERNAL_INCLUDE_DIR INTERFACE_INCLUDE AKANTU_INTERFACE_EXTERNAL_INCLUDE_DIR LIBRARIES AKANTU_EXTERNAL_LIBRARIES ) set(_include_dirs ${AKANTU_INCLUDE_DIRS} ${AKANTU_PRIVATE_EXTERNAL_INCLUDE_DIR} ${AKANTU_INTERFACE_EXTERNAL_INCLUDE_DIR} ) try_run(_material_list_run _material_list_compile ${CMAKE_BINARY_DIR} ${PROJECT_SOURCE_DIR}/cmake/material_lister.cc CMAKE_FLAGS "-DINCLUDE_DIRECTORIES:STRING=${_include_dirs}" "-DCMAKE_CXX_STANDARD=14" COMPILE_DEFINITIONS "-DAKANTU_CMAKE_LIST_MATERIALS" COMPILE_OUTPUT_VARIABLE _compile_results RUN_OUTPUT_VARIABLE _result_material_list) if(_material_list_compile AND "${_material_list_run}" EQUAL 0) message(STATUS "Materials included in Akantu:") string(REPLACE "\n" ";" _material_list "${_result_material_list}") foreach(_mat ${_material_list}) string(REPLACE ":" ";" _mat_key "${_mat}") list(GET _mat_key 0 _key) list(GET _mat_key 1 _class) list(LENGTH _mat_key _l) if("${_l}" GREATER 2) list(REMOVE_AT _mat_key 0 1) set(_opt " -- options: [") foreach(_o ${_mat_key}) set(_opt "${_opt} ${_o}") endforeach() set(_opt "${_opt} ]") else() set(_opt "") endif() message(STATUS " ${_class} -- key: ${_key}${_opt}") endforeach() else() message(STATUS "Could not determine the list of materials.") message("${_compile_results}") endif() endfunction() #=============================================================================== # Declare the options for the types and defines the approriate typedefs function(declare_akantu_types) set(AKANTU_TYPE_FLOAT "double (64bit)" CACHE STRING "Precision force floating point types") mark_as_advanced(AKANTU_TYPE_FLOAT) set_property(CACHE AKANTU_TYPE_FLOAT PROPERTY STRINGS "quadruple (128bit)" "double (64bit)" "float (32bit)" ) set(AKANTU_TYPE_INTEGER "int (32bit)" CACHE STRING "Size of the integer types") mark_as_advanced(AKANTU_TYPE_INTEGER) set_property(CACHE AKANTU_TYPE_INTEGER PROPERTY STRINGS "int (32bit)" "long int (64bit)" ) include(CheckTypeSize) # ---------------------------------------------------------------------------- # Floating point types # ---------------------------------------------------------------------------- if(AKANTU_TYPE_FLOAT STREQUAL "float (32bit)") set(AKANTU_FLOAT_TYPE "float" CACHE INTERNAL "") set(AKANTU_FLOAT_SIZE 4 CACHE INTERNAL "") elseif(AKANTU_TYPE_FLOAT STREQUAL "double (64bit)") set(AKANTU_FLOAT_TYPE "double" CACHE INTERNAL "") set(AKANTU_FLOAT_SIZE 8 CACHE INTERNAL "") elseif(AKANTU_TYPE_FLOAT STREQUAL "quadruple (128bit)") check_type_size("long double" LONG_DOUBLE) if(HAVE_LONG_DOUBLE) set(AKANTU_FLOAT_TYPE "long double" CACHE INTERNAL "") set(AKANTU_FLOAT_SIZE 16 CACHE INTERNAL "") message("This feature is not tested and will most probably not compile") else() message(FATAL_ERROR "The type long double is not defined on your system") endif() else() message(FATAL_ERROR "The float type is not defined") endif() include(CheckIncludeFileCXX) include(CheckCXXSourceCompiles) # ---------------------------------------------------------------------------- # Integer types # ---------------------------------------------------------------------------- check_include_file_cxx(cstdint HAVE_CSTDINT) if(NOT HAVE_CSTDINT) check_include_file_cxx(stdint.h HAVE_STDINT_H) if(HAVE_STDINT_H) list(APPEND _int_include stdint.h) endif() else() list(APPEND _int_include cstdint) endif() check_include_file_cxx(cstddef HAVE_CSTDDEF) if(NOT HAVE_CSTDINT) check_include_file_cxx(stddef.h HAVE_STDDEF_H) if(HAVE_STDINT_H) list(APPEND _int_include stddef.h) endif() else() list(APPEND _int_include cstddef) endif() if(AKANTU_TYPE_INTEGER STREQUAL "int (32bit)") set(AKANTU_INTEGER_SIZE 4 CACHE INTERNAL "") check_type_size("int" INT) if(INT EQUAL 4) set(AKANTU_SIGNED_INTEGER_TYPE "int" CACHE INTERNAL "") set(AKANTU_UNSIGNED_INTEGER_TYPE "unsigned int" CACHE INTERNAL "") else() check_type_size("int32_t" INT32_T LANGUAGE CXX) if(HAVE_INT32_T) set(AKANTU_SIGNED_INTEGER_TYPE "int32_t" CACHE INTERNAL "") set(AKANTU_UNSIGNED_INTEGER_TYPE "uint32_t" CACHE INTERNAL "") list(APPEND _extra_includes ${_int_include}) endif() endif() elseif(AKANTU_TYPE_INTEGER STREQUAL "long int (64bit)") set(AKANTU_INTEGER_SIZE 8 CACHE INTERNAL "") check_type_size("long int" LONG_INT) if(LONG_INT EQUAL 8) set(AKANTU_SIGNED_INTEGER_TYPE "long int" CACHE INTERNAL "") set(AKANTU_UNSIGNED_INTEGER_TYPE "unsigned long int" CACHE INTERNAL "") else() check_type_size("long long int" LONG_LONG_INT) if(HAVE_LONG_LONG_INT AND LONG_LONG_INT EQUAL 8) set(AKANTU_SIGNED_INTEGER_TYPE "long long int" CACHE INTERNAL "") set(AKANTU_UNSIGNED_INTEGER_TYPE "unsigned long long int" CACHE INTERNAL "") else() check_type_size("int64_t" INT64_T) if(HAVE_INT64_T) set(AKANTU_SIGNED_INTEGER_TYPE "int64_t" CACHE INTERNAL "") set(AKANTU_UNSIGNED_INTEGER_TYPE "uint64_t" CACHE INTERNAL "") list(APPEND _extra_includes ${_int_include}) endif() endif() endif() else() message(FATAL_ERROR "The integer type is not defined") endif() # ---------------------------------------------------------------------------- # includes # ---------------------------------------------------------------------------- foreach(_inc ${_extra_includes}) set(_incs "#include <${_inc}>\n${_incs}") endforeach() set(AKANTU_TYPES_EXTRA_INCLUDES ${_incs} CACHE INTERNAL "") endfunction() function(mask_package_options prefix) get_property(_list DIRECTORY PROPERTY VARIABLES) foreach(_var ${_list}) if("${_var}" MATCHES "^${prefix}.*") mark_as_advanced(${_var}) endif() endforeach() endfunction() diff --git a/src/model/common/time_step_solvers/time_step_solver_default.cc b/src/model/common/time_step_solvers/time_step_solver_default.cc index edf18bcc5..05950eaf6 100644 --- a/src/model/common/time_step_solvers/time_step_solver_default.cc +++ b/src/model/common/time_step_solvers/time_step_solver_default.cc @@ -1,334 +1,332 @@ /** * @file time_step_solver_default.cc * * @author Nicolas Richart * * @date creation: Tue Sep 15 2015 * @date last modification: Wed Feb 21 2018 * * @brief Default implementation of the time step solver * * @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 "time_step_solver_default.hh" #include "dof_manager_default.hh" #include "integration_scheme_1st_order.hh" #include "integration_scheme_2nd_order.hh" #include "mesh.hh" #include "non_linear_solver.hh" #include "pseudo_time.hh" #include "sparse_matrix_aij.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ TimeStepSolverDefault::TimeStepSolverDefault( DOFManager & dof_manager, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver, SolverCallback & solver_callback, const ID & id, UInt memory_id) : TimeStepSolver(dof_manager, type, non_linear_solver, solver_callback, id, memory_id) { switch (type) { case TimeStepSolverType::_dynamic: break; case TimeStepSolverType::_dynamic_lumped: this->is_mass_lumped = true; break; case TimeStepSolverType::_static: /// initialize a static time solver for callback dofs break; default: AKANTU_TO_IMPLEMENT(); } } /* -------------------------------------------------------------------------- */ void TimeStepSolverDefault::setIntegrationSchemeInternal( const ID & dof_id, const IntegrationSchemeType & type, IntegrationScheme::SolutionType solution_type) { if (this->integration_schemes.find(dof_id) != this->integration_schemes.end()) { AKANTU_EXCEPTION("Their DOFs " << dof_id << " have already an integration scheme associated"); } std::unique_ptr integration_scheme; if (this->is_mass_lumped) { switch (type) { case IntegrationSchemeType::_forward_euler: { integration_scheme = std::make_unique(_dof_manager, dof_id); break; } case IntegrationSchemeType::_central_difference: { integration_scheme = std::make_unique(_dof_manager, dof_id); break; } default: AKANTU_EXCEPTION( "This integration scheme cannot be used in lumped dynamic"); } } else { switch (type) { case IntegrationSchemeType::_pseudo_time: { integration_scheme = std::make_unique(_dof_manager, dof_id); break; } case IntegrationSchemeType::_forward_euler: { integration_scheme = std::make_unique(_dof_manager, dof_id); break; } case IntegrationSchemeType::_trapezoidal_rule_1: { integration_scheme = std::make_unique(_dof_manager, dof_id); break; } case IntegrationSchemeType::_backward_euler: { integration_scheme = std::make_unique(_dof_manager, dof_id); break; } case IntegrationSchemeType::_central_difference: { integration_scheme = std::make_unique(_dof_manager, dof_id); break; } case IntegrationSchemeType::_fox_goodwin: { integration_scheme = std::make_unique(_dof_manager, dof_id); break; } case IntegrationSchemeType::_trapezoidal_rule_2: { integration_scheme = std::make_unique(_dof_manager, dof_id); break; } case IntegrationSchemeType::_linear_acceleration: { integration_scheme = std::make_unique(_dof_manager, dof_id); break; } case IntegrationSchemeType::_generalized_trapezoidal: { integration_scheme = std::make_unique(_dof_manager, dof_id); break; } case IntegrationSchemeType::_newmark_beta: integration_scheme = std::make_unique(_dof_manager, dof_id); break; } } AKANTU_DEBUG_ASSERT(integration_scheme, "No integration scheme was found for the provided types"); auto && matrices_names = integration_scheme->getNeededMatrixList(); for (auto && name : matrices_names) { needed_matrices.insert({name, _mt_not_defined}); } this->integration_schemes[dof_id] = std::move(integration_scheme); this->solution_types[dof_id] = solution_type; this->integration_schemes_owner.insert(dof_id); } /* -------------------------------------------------------------------------- */ bool TimeStepSolverDefault::hasIntegrationScheme(const ID & dof_id) const { return this->integration_schemes.find(dof_id) != this->integration_schemes.end(); } /* -------------------------------------------------------------------------- */ TimeStepSolverDefault::~TimeStepSolverDefault() = default; /* -------------------------------------------------------------------------- */ void TimeStepSolverDefault::solveStep(SolverCallback & solver_callback) { this->solver_callback = &solver_callback; - this->non_linear_solver.solve(*this); - this->solver_callback = nullptr; } /* -------------------------------------------------------------------------- */ void TimeStepSolverDefault::predictor() { TimeStepSolver::predictor(); for (auto && pair : this->integration_schemes) { const auto & dof_id = pair.first; auto & integration_scheme = pair.second; if (this->_dof_manager.hasPreviousDOFs(dof_id)) { this->_dof_manager.savePreviousDOFs(dof_id); } /// integrator predictor integration_scheme->predictor(this->time_step); } } /* -------------------------------------------------------------------------- */ void TimeStepSolverDefault::corrector() { AKANTU_DEBUG_IN(); TimeStepSolver::corrector(); for (auto & pair : this->integration_schemes) { auto & dof_id = pair.first; auto & integration_scheme = pair.second; const auto & solution_type = this->solution_types[dof_id]; integration_scheme->corrector(solution_type, this->time_step); /// computing the increment of dof if needed if (this->_dof_manager.hasDOFsIncrement(dof_id)) { if (not this->_dof_manager.hasPreviousDOFs(dof_id)) { AKANTU_DEBUG_WARNING("In order to compute the increment of " << dof_id << " a 'previous' has to be registered"); continue; } auto & increment = this->_dof_manager.getDOFsIncrement(dof_id); auto & previous = this->_dof_manager.getPreviousDOFs(dof_id); auto dof_array_comp = this->_dof_manager.getDOFs(dof_id).getNbComponent(); increment.copy(this->_dof_manager.getDOFs(dof_id)); for (auto && data : zip(make_view(increment, dof_array_comp), make_view(previous, dof_array_comp))) { std::get<0>(data) -= std::get<1>(data); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void TimeStepSolverDefault::assembleMatrix(const ID & matrix_id) { AKANTU_DEBUG_IN(); TimeStepSolver::assembleMatrix(matrix_id); if (matrix_id != "J") return; for_each_integrator([&](auto && dof_id, auto && integration_scheme) { const auto & solution_type = this->solution_types[dof_id]; integration_scheme.assembleJacobian(solution_type, this->time_step); }); this->_dof_manager.applyBoundary("J"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ // void TimeStepSolverDefault::assembleLumpedMatrix(const ID & matrix_id) { // AKANTU_DEBUG_IN(); // TimeStepSolver::assembleLumpedMatrix(matrix_id); // if (matrix_id != "J") // return; // for (auto & pair : this->integration_schemes) { // auto & dof_id = pair.first; // auto & integration_scheme = pair.second; // const auto & solution_type = this->solution_types[dof_id]; // integration_scheme->assembleJacobianLumped(solution_type, // this->time_step); // } // this->_dof_manager.applyBoundaryLumped("J"); // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ void TimeStepSolverDefault::assembleResidual() { if (this->needed_matrices.find("M") != needed_matrices.end()) { if (this->is_mass_lumped) { this->assembleLumpedMatrix("M"); } else { this->assembleMatrix("M"); } } TimeStepSolver::assembleResidual(); for_each_integrator([&](auto &&, auto && integration_scheme) { integration_scheme.assembleResidual(this->is_mass_lumped); }); } /* -------------------------------------------------------------------------- */ void TimeStepSolverDefault::assembleResidual(const ID & residual_part) { AKANTU_DEBUG_IN(); if (this->needed_matrices.find("M") != needed_matrices.end()) { if (this->is_mass_lumped) { this->assembleLumpedMatrix("M"); } else { this->assembleMatrix("M"); } } if (residual_part != "inertial") { TimeStepSolver::assembleResidual(residual_part); } if (residual_part == "inertial") { for_each_integrator([&](auto &&, auto && integration_scheme) { integration_scheme.assembleResidual(this->is_mass_lumped); }); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void TimeStepSolverDefault::beforeSolveStep() { TimeStepSolver::beforeSolveStep(); for_each_integrator( [&](auto &&, auto && integration_scheme) { integration_scheme.store(); }); } /* -------------------------------------------------------------------------- */ void TimeStepSolverDefault::afterSolveStep(bool converged) { if (not converged) { for_each_integrator([&](auto &&, auto && integration_scheme) { integration_scheme.restore(); }); } TimeStepSolver::afterSolveStep(converged); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc index 50b15e2f9..c568c45f6 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model.cc @@ -1,1200 +1,1201 @@ /** * @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), material_index("material index", id, memory_id), material_local_numbering("material local numbering", id, memory_id) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("SolidMechanicsFEEngine", mesh, Model::spatial_dimension); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("solid_mechanics_model", 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() = default; /* -------------------------------------------------------------------------- */ 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->assembleInternalForces(); 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; if (matrix_id == "K") { auto matrix_type = _unsymmetric; for (auto & material : materials) { matrix_type = std::max(matrix_type, material->getMatrixType(matrix_id)); } } 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(bool converged) { for (auto & material : materials) material->afterSolveStep(converged); } /* -------------------------------------------------------------------------- */ 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->hasMatrixChanged("K"); } 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 (mesh.getSpatialDimension(elem.type) != spatial_dimension) continue; 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 fails in parallel if the event is sent on facet between constructor // and initFull \todo: to debug... 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 (not aka::is_of_type(*mat)) 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 (not aka::is_of_type(*mat)) 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 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/test/test_model/patch_tests/patch_test_linear_elastic_explicit.cc b/test/test_model/patch_tests/patch_test_linear_elastic_explicit.cc index f5872b3a9..851d82459 100644 --- a/test/test_model/patch_tests/patch_test_linear_elastic_explicit.cc +++ b/test/test_model/patch_tests/patch_test_linear_elastic_explicit.cc @@ -1,85 +1,98 @@ /** * @file patch_test_linear_elastic_explicit.cc * * @author Nicolas Richart * * @date creation: Tue Jan 30 2018 * * @brief patch test solid mechanics explicit * * @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 "patch_test_linear_solid_mechanics_fixture.hh" /* -------------------------------------------------------------------------- */ TYPED_TEST(TestPatchTestSMMLinear, Explicit) { std::string filename = "material_check_stress_plane_stress.dat"; if (this->plane_strain) filename = "material_check_stress_plane_strain.dat"; this->initModel(_explicit_lumped_mass, filename); const auto & coordinates = this->mesh->getNodes(); auto & displacement = this->model->getDisplacement(); // set the position of all nodes to the static solution for (auto && tuple : zip(make_view(coordinates, this->dim), make_view(displacement, this->dim))) { this->setLinearDOF(std::get<1>(tuple), std::get<0>(tuple)); } for (UInt s = 0; s < 100; ++s) { this->model->solveStep(); } auto ekin = this->model->getEnergy("kinetic"); EXPECT_NEAR(0, ekin, 1e-16); this->checkAll(); + +#define debug 0 +#if debug + this->model->setBaseName(std::to_string(this->type) + "_explicit"); + this->model->addDumpField("stress"); + this->model->addDumpField("grad_u"); + this->model->addDumpFieldVector("internal_force"); + this->model->addDumpFieldVector("external_force"); + this->model->addDumpField("blocked_dofs"); + this->model->addDumpFieldVector("displacement"); + this->model->dump(); +#endif + } /* -------------------------------------------------------------------------- */ TYPED_TEST(TestPatchTestSMMLinear, ExplicitFiniteDeformation) { std::string filename = "material_check_stress_plane_stress_finite_deformation.dat"; if (this->plane_strain) filename = "material_check_stress_plane_strain_finite_deformation.dat"; else { SUCCEED(); return; } this->initModel(_explicit_lumped_mass, filename); const auto & coordinates = this->mesh->getNodes(); auto & displacement = this->model->getDisplacement(); // set the position of all nodes to the static solution for (auto && tuple : zip(make_view(coordinates, this->dim), make_view(displacement, this->dim))) { this->setLinearDOF(std::get<1>(tuple), std::get<0>(tuple)); } for (UInt s = 0; s < 100; ++s) { this->model->solveStep(); } auto ekin = this->model->getEnergy("kinetic"); EXPECT_NEAR(0, ekin, 1e-16); this->checkAll(); } diff --git a/test/test_model/patch_tests/patch_test_linear_elastic_implicit.cc b/test/test_model/patch_tests/patch_test_linear_elastic_implicit.cc index 3e55bfa9e..c18145556 100644 --- a/test/test_model/patch_tests/patch_test_linear_elastic_implicit.cc +++ b/test/test_model/patch_tests/patch_test_linear_elastic_implicit.cc @@ -1,130 +1,170 @@ /** * @file patch_test_linear_elastic_implicit.cc * * @author Nicolas Richart * * @date creation: Tue Jan 30 2018 * * @brief Patch test for SolidMechanics implicit * * @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 "patch_test_linear_solid_mechanics_fixture.hh" /* -------------------------------------------------------------------------- */ #include "non_linear_solver.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; TYPED_TEST(TestPatchTestSMMLinear, Implicit) { std::string filename = "material_check_stress_plane_stress.dat"; if (this->plane_strain) filename = "material_check_stress_plane_strain.dat"; this->initModel(_implicit_dynamic, filename); const auto & coordinates = this->mesh->getNodes(); auto & displacement = this->model->getDisplacement(); // set the position of all nodes to the static solution for (auto && tuple : zip(make_view(coordinates, this->dim), make_view(displacement, this->dim))) { this->setLinearDOF(std::get<1>(tuple), std::get<0>(tuple)); } for (UInt s = 0; s < 100; ++s) { this->model->solveStep(); } auto ekin = this->model->getEnergy("kinetic"); EXPECT_NEAR(0, ekin, 1e-16); this->checkAll(); + +#define debug 0 +#if debug + this->model->setBaseName(std::to_string(this->type) + "_implicit"); + this->model->addDumpField("stress"); + this->model->addDumpField("grad_u"); + this->model->addDumpFieldVector("internal_force"); + this->model->addDumpFieldVector("external_force"); + this->model->addDumpField("blocked_dofs"); + this->model->addDumpFieldVector("displacement"); + this->model->dump(); +#endif + } /* -------------------------------------------------------------------------- */ TYPED_TEST(TestPatchTestSMMLinear, Static) { std::string filename = "material_check_stress_plane_stress.dat"; if (this->plane_strain) filename = "material_check_stress_plane_strain.dat"; this->initModel(_static, filename); auto & solver = this->model->getNonLinearSolver(); solver.set("max_iterations", 2); solver.set("threshold", 2e-4); solver.set("convergence_type", SolveConvergenceCriteria::_residual); this->model->solveStep(); this->checkAll(); + +#define debug 0 +#if debug + this->model->setBaseName(std::to_string(this->type) + "_static"); + this->model->addDumpField("stress"); + this->model->addDumpField("grad_u"); + this->model->addDumpFieldVector("internal_force"); + this->model->addDumpFieldVector("external_force"); + this->model->addDumpField("blocked_dofs"); + this->model->addDumpFieldVector("displacement"); + this->model->dump(); +#endif + } /* -------------------------------------------------------------------------- */ TYPED_TEST(TestPatchTestSMMLinear, ImplicitFiniteDeformation) { std::string filename = "material_check_stress_plane_stress_finite_deformation.dat"; if (this->plane_strain) filename = "material_check_stress_plane_strain_finite_deformation.dat"; else { SUCCEED(); return; } this->initModel(_implicit_dynamic, filename); const auto & coordinates = this->mesh->getNodes(); auto & displacement = this->model->getDisplacement(); // set the position of all nodes to the static solution for (auto && tuple : zip(make_view(coordinates, this->dim), make_view(displacement, this->dim))) { this->setLinearDOF(std::get<1>(tuple), std::get<0>(tuple)); } for (UInt s = 0; s < 100; ++s) { this->model->solveStep(); } auto ekin = this->model->getEnergy("kinetic"); EXPECT_NEAR(0, ekin, 1e-16); this->checkAll(); + + +#define debug 0 +#if debug + this->model->setBaseName(std::to_string(this->type) + "_implicit_finit_def"); + this->model->addDumpField("stress"); + this->model->addDumpField("grad_u"); + this->model->addDumpFieldVector("internal_force"); + this->model->addDumpFieldVector("external_force"); + this->model->addDumpField("blocked_dofs"); + this->model->addDumpFieldVector("displacement"); + this->model->dump(); +#endif + } /* -------------------------------------------------------------------------- */ TYPED_TEST(TestPatchTestSMMLinear, StaticFiniteDeformation) { std::string filename = "material_check_stress_plane_stress_finite_deformation.dat"; if (this->plane_strain) { filename = "material_check_stress_plane_strain_finite_deformation.dat"; } else { SUCCEED(); return; } this->initModel(_static, filename); auto & solver = this->model->getNonLinearSolver(); solver.set("max_iterations", 2); solver.set("threshold", 2e-4); solver.set("convergence_type", SolveConvergenceCriteria::_residual); this->model->solveStep(); this->checkAll(); } diff --git a/test/test_model/patch_tests/patch_test_linear_solid_mechanics_fixture.hh b/test/test_model/patch_tests/patch_test_linear_solid_mechanics_fixture.hh index b92871ceb..fb5b5cce0 100644 --- a/test/test_model/patch_tests/patch_test_linear_solid_mechanics_fixture.hh +++ b/test/test_model/patch_tests/patch_test_linear_solid_mechanics_fixture.hh @@ -1,113 +1,154 @@ /** * @file patch_test_linear_solid_mechanics_fixture.hh * * @author Nicolas Richart * * @date creation: Tue Jan 30 2018 * * @brief SolidMechanics patch tests fixture * * @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 "patch_test_linear_fixture.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_PATCH_TEST_LINEAR_SOLID_MECHANICS_FIXTURE_HH__ #define __AKANTU_PATCH_TEST_LINEAR_SOLID_MECHANICS_FIXTURE_HH__ /* -------------------------------------------------------------------------- */ template class TestPatchTestSMMLinear : public TestPatchTestLinear, SolidMechanicsModel> { using parent = TestPatchTestLinear, SolidMechanicsModel>; public: static constexpr bool plane_strain = std::tuple_element_t<1, tuple_>::value; void applyBC() override { parent::applyBC(); auto & displacement = this->model->getDisplacement(); this->applyBConDOFs(displacement); } + void checkForces() { + auto & mat = this->model->getMaterial(0); + auto & internal_forces = this->model->getInternalForce(); + auto & external_forces = this->model->getExternalForce(); + auto dim = this->dim; + + Matrix sigma = + make_view(mat.getStress(this->type), dim, dim).begin()[0]; + + external_forces.clear(); + if (dim > 1) { + for (auto & eg : this->mesh->iterateElementGroups()) { + this->model->applyBC(BC::Neumann::FromHigherDim(sigma), eg.getName()); + } + } else { + external_forces(0) = -sigma(0, 0); + external_forces(1) = sigma(0, 0); + } + + Real force_norm_inf = -std::numeric_limits::max(); + + Vector total_force(dim); + total_force.clear(); + + for (auto && f : make_view(internal_forces, dim)) { + total_force += f; + force_norm_inf = std::max(force_norm_inf, f.template norm()); + } + + EXPECT_NEAR(0, total_force.template norm() / force_norm_inf, 1e-9); + + for (auto && tuple : zip(make_view(internal_forces, dim), + make_view(external_forces, dim))) { + auto && f_int = std::get<0>(tuple); + auto && f_ext = std::get<1>(tuple); + auto f = f_int + f_ext; + EXPECT_NEAR(0, f.template norm() / force_norm_inf, 1e-9); + } + } + void checkAll() { auto & displacement = this->model->getDisplacement(); auto & mat = this->model->getMaterial(0); this->checkDOFs(displacement); this->checkGradient(mat.getGradU(this->type), displacement); this->checkResults( [&](const Matrix & pstrain) { Real nu = this->model->getMaterial(0).get("nu"); Real E = this->model->getMaterial(0).get("E"); auto strain = (pstrain + pstrain.transpose()) / 2.; auto trace = strain.trace(); auto lambda = nu * E / ((1 + nu) * (1 - 2 * nu)); auto mu = E / (2 * (1 + nu)); if (not this->plane_strain) { lambda = nu * E / (1 - nu * nu); } decltype(strain) stress(this->dim, this->dim); if (this->dim == 1) { stress(0, 0) = E * strain(0, 0); } else { for (UInt i = 0; i < this->dim; ++i) for (UInt j = 0; j < this->dim; ++j) stress(i, j) = (i == j) * lambda * trace + 2 * mu * strain(i, j); } return stress; }, mat.getStress(this->type), displacement); + this->checkForces(); } }; template constexpr bool TestPatchTestSMMLinear::plane_strain; template struct invalid_plan_stress : std::true_type {}; template struct invalid_plan_stress> : aka::bool_constant::getSpatialDimension() != 2 and not bool_c::value> {}; using true_false = std::tuple, aka::bool_constant>; template using valid_types = aka::negation>; using model_types = gtest_list_t< tuple_filter_t>>; TYPED_TEST_SUITE(TestPatchTestSMMLinear, model_types); #endif /* __AKANTU_PATCH_TEST_LINEAR_SOLID_MECHANICS_FIXTURE_HH__ */