diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 65747616a..ab607bb50 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -1,350 +1,345 @@ # yaml-language-server: $schema=gitlab-ci # yaml-language-server: $format.enable=false stages: - configure - build_libs - build_tests - test - code_quality - deploy include: local: '.gitlab-ci.d/*.yaml' #------------------------------------------------------------------------------- # Rebuilding the docker images if needed #------------------------------------------------------------------------------- docker build:debian-bullseye: variables: IMAGE_NAME: debian:bullseye extends: .docker_build docker build:ubuntu-lts: variables: IMAGE_NAME: ubuntu:lts extends: .docker_build docker build:manylinux: variables: IMAGE_NAME: manylinux:2010_x86_64 extends: .docker_build # ------------------------------------------------------------------------------ # Debian bullseye compiled with GCC # ------------------------------------------------------------------------------ configure:debian_bullseye_gcc: extends: - .debian_bullseye_gcc - .build_coverage - .configure build:debian_bullseye_gcc: extends: - .debian_bullseye_gcc - .build_coverage - .build_all needs: - job: configure:debian_bullseye_gcc test:debian_bullseye_gcc: extends: - .debian_bullseye_gcc - .build_coverage - .tests coverage: '/^lines: (\d+\.\d+\%)/' needs: - job: build:debian_bullseye_gcc # ------------------------------------------------------------------------------ # Debian bullseye compiled with Clang # ------------------------------------------------------------------------------ configure:debian_bullseye_clang: extends: - .debian_bullseye_clang - .build_coverage - .configure build:debian_bullseye_clang: extends: - .debian_bullseye_clang - .build_coverage - .build_all needs: - job: configure:debian_bullseye_clang test:debian_bullseye_clang: extends: - .debian_bullseye_clang - .build_coverage - .tests coverage: '/^lines: (\d+\.\d+\%)/' needs: - job: build:debian_bullseye_clang # ------------------------------------------------------------------------------ # Ubuntu LTS compiled with GCC # ------------------------------------------------------------------------------ configure:ubuntu_lts_gcc: extends: - .ubuntu_lts_gcc - .build_release - .configure build:ubuntu_lts_gcc: extends: - .ubuntu_lts_gcc - .build_release - .build_all needs: - job: configure:ubuntu_lts_gcc test:ubuntu_lts_gcc: extends: - .ubuntu_lts_gcc - .build_release - .tests needs: - job: build:ubuntu_lts_gcc # ------------------------------------------------------------------------------ # Debian bullseye compiled with GCC tested with valgrind # ------------------------------------------------------------------------------ configure:ubuntu_lts_gcc_valgrind: extends: - .ubuntu_lts_gcc - .build_valgrind - .configure build:ubuntu_lts_gcc_valgrind: extends: - .ubuntu_lts_gcc - .build_valgrind - .build_all needs: - job: configure:ubuntu_lts_gcc_valgrind test:ubuntu_lts_gcc_valgrind: extends: - .ubuntu_lts_gcc - .build_valgrind - .tests needs: - job: build:ubuntu_lts_gcc_valgrind # ------------------------------------------------------------------------------ # Manylinux to build python packages # ------------------------------------------------------------------------------ configure:python_package: stage: configure extends: - .manylinux_2010_x64_gcc - .build_release script: # create the build folder - cmake -E make_directory build - cd build # Variables for cmake - export CMAKE_PREFIX_PATH=/softs/view - export BOOST_ROOT=/softs/view # Configure in sequential and without tests or examples - cmake -DAKANTU_COHESIVE_ELEMENT:BOOL=TRUE -DAKANTU_IMPLICIT:BOOL=TRUE -DAKANTU_PARALLEL:BOOL=FALSE -DAKANTU_STRUCTURAL_MECHANICS:BOOL=TRUE -DAKANTU_HEAT_TRANSFER:BOOL=TRUE -DAKANTU_DAMAGE_NON_LOCAL:BOOL=TRUE -DAKANTU_PHASE_FIELD:BOOL=TRUE -DAKANTU_PYTHON_INTERFACE:BOOL=FALSE -DAKANTU_CONTACT_MECHANICS:BOOL=TRUE -DAKANTU_EXAMPLES:BOOL=FALSE -DAKANTU_TESTS:BOOL=FALSE -DMUMPS_DETECT_DEBUG:BOOL=TRUE -DCMAKE_INSTALL_PREFIX:PATH=${CI_PROJECT_DIR}/install -DCMAKE_BUILD_TYPE:STRING=${BUILD_TYPE} .. artifacts: when: on_success paths: - build/ expire_in: 10h build_akantu:python_package: extends: - .build_libs - .build_release - .manylinux_2010_x64_gcc - script: - stage: build_libs script: - cmake --build build --target akantu -j1 - cmake --install build artifacts: when: on_success paths: - install/ expire_in: 10h needs: - job: configure:python_package build_pip:python_package: stage: build_tests extends: - .build_release - .manylinux_2010_x64_gcc script: - export CI_AKANTU_INSTALL_PREFIX=${CI_PROJECT_DIR}/install - export CMAKE_PREFIX_PATH=/softs/view:${CI_AKANTU_INSTALL_PREFIX} - test/ci/make-wheels.sh needs: - job: build_akantu:python_package artifacts: when: on_success paths: - wheelhouse expire_in: 10h test:python_package: stage: test image: python:3.8 needs: - job: build_pip:python_package script: - pip install numpy scipy - pip install akantu --no-index --find-links=${PWD}/wheelhouse - python -c "import akantu" - cd examples/python/dynamics/ - apt update && apt install -y gmsh - gmsh -2 bar.geo - python ./dynamics.py package:python_gitlab: stage: deploy image: python:latest script: - pip install twine - TWINE_PASSWORD=${CI_JOB_TOKEN} TWINE_USERNAME=gitlab-ci-token python3 -m twine upload --repository-url https://gitlab.com/api/v4/projects/${CI_PROJECT_ID}/packages/pypi wheelhouse/* needs: - job: build_pip:python_package - job: test:python_package only: - master package:python_pypi: stage: deploy image: python:latest script: - pip install twine - TWINE_PASSWORD=${PYPI_TOKEN} TWINE_USERNAME=__token__ python3 -m twine upload --verbose wheelhouse/* needs: - job: build_pip:python_package - job: test:python_package only: - tags # ------------------------------------------------------------------------------ # Code Quality # ------------------------------------------------------------------------------ cq:code_quality: extends: - .code_quality_gitlab_template needs: - job: build:debian_bullseye_clang artifacts: paths: - gl-code-quality-report.json cq:clang_tidy: extends: - .clang_tools script: - test/ci/scripts/cq -x third-party -x extra-packages -x pybind11 -x test ${FILE_LIST_ARG} clang-tidy -p ${CI_PROJECT_DIR}/build > gl-clang-tidy-report.json needs: - job: build:debian_bullseye_clang artifacts: paths: - gl-clang-tidy-report.json cq:clang_format: extends: - .clang_tools script: - test/ci/scripts/cq -x third-party -x extra-packages clang-format -p ${CI_PROJECT_DIR}/build > gl-clang-format-report.json needs: - job: build:debian_bullseye_clang artifacts: paths: - gl-clang-format-report.json cq:compilation_warnings: stage: code_quality image: python:latest script: - pip install warning-parser termcolor Click - ls build-*-err.log - test/ci/scripts/cq -x third-party -x extra-packages warnings build-*-err.log > gl-warnings-report.json needs: - job: build:debian_bullseye_clang - job: build:debian_bullseye_gcc - job: build:ubuntu_lts_gcc artifacts: paths: - gl-warnings-report.json cq:merge_code_quality: stage: deploy extends: - .debian_bullseye_clang script: - jq -Ms '[.[][]]' gl-*-report.json | tee gl-codequality.json | jq -C needs: - job: cq:code_quality - job: cq:clang_tidy - job: cq:clang_format - job: cq:compilation_warnings - artifacts: - paths: - - gl-codequality.json artifacts: reports: codequality: [gl-codequality.json] # ------------------------------------------------------------------------------ # Deploy pages # ------------------------------------------------------------------------------ pages: stage: deploy extends: - .debian_bullseye_gcc script: - cd build - cmake -DAKANTU_DOCUMENTATION=ON .. - cmake --build . -t sphinx-doc - mv doc/dev-doc/html ../public needs: - job: build:debian_bullseye_gcc artifacts: paths: - public only: - master diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc index 21dc28321..f33020eb4 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model.cc @@ -1,1250 +1,1254 @@ /** * @file solid_mechanics_model.cc * * @author Ramin Aghababaei * @author Guillaume Anciaux * @author Mauro Corrado * @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: Fri Apr 09 2021 * * @brief Implementation of the SolidMechanicsModel class * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "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" /* -------------------------------------------------------------------------- */ #include "dumper_iohelper_paraview.hh" /* -------------------------------------------------------------------------- */ #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 * @param model_type this is an internal parameter for inheritance purposes */ SolidMechanicsModel::SolidMechanicsModel( Mesh & mesh, UInt dim, const ID & id, std::shared_ptr dof_manager, const ModelType model_type) : Model(mesh, model_type, std::move(dof_manager), dim, id), material_index("material index", id), material_local_numbering("material local numbering", id) { AKANTU_DEBUG_IN(); this->registerFEEngineObject("SolidMechanicsFEEngine", mesh, Model::spatial_dimension); this->mesh.registerDumper("solid_mechanics_model", id, true); this->mesh.addDumpMesh(mesh, Model::spatial_dimension, _not_ghost, _ek_regular); material_selector = std::make_shared(material_index); 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); this->mesh.getDumper().setTimeStep(time_step); } /* -------------------------------------------------------------------------- */ /* 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 options * \parblock * contains the different options to initialize the model * \li \c analysis_method specify the type of solver to use * \endparblock */ 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 (not this->parser.getLastParsedFile().empty()) { 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 /*unused*/) { 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) const { // \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$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->zero(); // 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(bool need_to_reassemble) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix."); // Check if materials need to recompute the matrix for (auto & material : materials) { need_to_reassemble |= material->hasMatrixChanged("K"); } if (need_to_reassemble) { this->getDOFManager().getMatrix("K").zero(); // 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(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")) { + this->assembleLumpedMatrix("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")) { + this->assembleMatrix("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)) * static_cast(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(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(); } 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, 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; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const ID & energy_id, const ID & group_id) { auto && group = mesh.getElementGroup(group_id); auto energy = 0.; for (auto && type : group.elementTypes()) { for (auto el : group.getElementsIterable(type)) { energy += getEnergy(energy_id, el); } } /// reduction sum over all processors mesh.getCommunicator().allReduce(energy, SynchronizerOperation::_sum); 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()); for (const 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, false); } 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 (const 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( 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); quadrature_points_coordinates.initialize(this->getFEEngine(), _nb_component = spatial_dimension, _ghost_type = ghost_type); for (const 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(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, GhostType ghost_type, 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, GhostType ghost_type, 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: { 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: /* FALLTHRU */ 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/structural_mechanics/structural_mechanics_model.cc b/src/model/structural_mechanics/structural_mechanics_model.cc index 7e111852b..6c01e0cb1 100644 --- a/src/model/structural_mechanics/structural_mechanics_model.cc +++ b/src/model/structural_mechanics/structural_mechanics_model.cc @@ -1,786 +1,777 @@ /** * @file structural_mechanics_model.cc * * @author Fabian Barras * @author Lucas Frerot * @author Sébastien Hartmann * @author Nicolas Richart * @author Damien Spielmann * * @date creation: Fri Jul 15 2011 * @date last modification: Mon Mar 15 2021 * * @brief Model implementation for Structural Mechanics elements * * * @section LICENSE * * Copyright (©) 2010-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "structural_mechanics_model.hh" #include "dof_manager.hh" #include "integrator_gauss.hh" #include "mesh.hh" #include "shape_structural.hh" #include "sparse_matrix.hh" #include "time_step_solver.hh" /* -------------------------------------------------------------------------- */ #include "dumpable_inline_impl.hh" #include "dumper_elemental_field.hh" #include "dumper_internal_material_field.hh" #include "dumper_iohelper_paraview.hh" #include "group_manager_inline_impl.hh" /* -------------------------------------------------------------------------- */ #include "structural_element_bernoulli_beam_2.hh" #include "structural_element_bernoulli_beam_3.hh" #include "structural_element_kirchhoff_shell.hh" /* -------------------------------------------------------------------------- */ //#include "structural_mechanics_model_inline_impl.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ inline UInt StructuralMechanicsModel::getNbDegreeOfFreedom(ElementType type) { UInt ndof = 0; #define GET_(type) ndof = ElementClass::getNbDegreeOfFreedom() AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_, _ek_structural); #undef GET_ return ndof; } /* -------------------------------------------------------------------------- */ StructuralMechanicsModel::StructuralMechanicsModel(Mesh & mesh, UInt dim, const ID & id) : Model(mesh, ModelType::_structural_mechanics_model, dim, id), f_m2a(1.0), stress("stress", id), element_material("element_material", id), set_ID("beam sets", id) { AKANTU_DEBUG_IN(); registerFEEngineObject("StructuralMechanicsFEEngine", mesh, spatial_dimension); if (spatial_dimension == 2) { nb_degree_of_freedom = 3; } else if (spatial_dimension == 3) { nb_degree_of_freedom = 6; } else { AKANTU_TO_IMPLEMENT(); } this->mesh.registerDumper("structural_mechanics_model", id, true); this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_structural); this->initDOFManager(); this->dumper_default_element_kind = _ek_structural; mesh.getElementalData("extra_normal") .initialize(mesh, _element_kind = _ek_structural, _nb_component = spatial_dimension, _with_nb_element = true, _default_value = 0.); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ StructuralMechanicsModel::~StructuralMechanicsModel() = default; /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initFullImpl(const ModelOptions & options) { Model::initFullImpl(options); // Initializing stresses ElementTypeMap stress_components; /// TODO this is ugly af, maybe add a function to FEEngine for (auto && type : mesh.elementTypes(_spatial_dimension = _all_dimensions, _element_kind = _ek_structural)) { UInt nb_components = 0; // Getting number of components for each element type #define GET_(type) nb_components = ElementClass::getNbStressComponents() AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(GET_); #undef GET_ stress_components(nb_components, type); } stress.initialize( getFEEngine(), _spatial_dimension = _all_dimensions, _element_kind = _ek_structural, _nb_component = [&stress_components](ElementType type, GhostType /*unused*/) -> UInt { return stress_components(type); }); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initFEEngineBoundary() { /// TODO: this function should not be reimplemented /// we're just avoiding a call to Model::initFEEngineBoundary() } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::setTimeStep(Real time_step, const ID & solver_id) { Model::setTimeStep(time_step, solver_id); this->mesh.getDumper().setTimeStep(time_step); } /* -------------------------------------------------------------------------- */ /* Initialisation */ /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initSolver( TimeStepSolverType time_step_solver_type, NonLinearSolverType /*unused*/) { AKANTU_DEBUG_IN(); this->allocNodalField(displacement_rotation, nb_degree_of_freedom, "displacement"); this->allocNodalField(external_force, nb_degree_of_freedom, "external_force"); this->allocNodalField(internal_force, nb_degree_of_freedom, "internal_force"); this->allocNodalField(blocked_dofs, nb_degree_of_freedom, "blocked_dofs"); auto & dof_manager = this->getDOFManager(); if (not dof_manager.hasDOFs("displacement")) { dof_manager.registerDOFs("displacement", *displacement_rotation, _dst_nodal); dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs); } if (time_step_solver_type == TimeStepSolverType::_dynamic || time_step_solver_type == TimeStepSolverType::_dynamic_lumped) { this->allocNodalField(velocity, nb_degree_of_freedom, "velocity"); this->allocNodalField(acceleration, nb_degree_of_freedom, "acceleration"); if (!dof_manager.hasDOFsDerivatives("displacement", 1)) { dof_manager.registerDOFsDerivative("displacement", 1, *this->velocity); dof_manager.registerDOFsDerivative("displacement", 2, *this->acceleration); } /* Only allocate the mass if the "lumped" mode is ennabled. * Also it is not a 1D array, but has an element for every * DOF, which are most of the time equal, but makes handling * some operations a bit simpler. */ if (time_step_solver_type == TimeStepSolverType::_dynamic_lumped) { this->allocateLumpedMassArray(); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initModel() { element_material.initialize(mesh, _element_kind = _ek_structural, _default_value = 0, _with_nb_element = true); getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); if (not need_to_reassemble_stiffness) { return; } if (not getDOFManager().hasMatrix("K")) { getDOFManager().getNewMatrix("K", getMatrixType("K")); } this->getDOFManager().zeroMatrix("K"); for (const auto & type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_structural)) { #define ASSEMBLE_STIFFNESS_MATRIX(type) assembleStiffnessMatrix(); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(ASSEMBLE_STIFFNESS_MATRIX); #undef ASSEMBLE_STIFFNESS_MATRIX } need_to_reassemble_stiffness = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::computeStresses() { AKANTU_DEBUG_IN(); for (const auto & type : mesh.elementTypes(spatial_dimension, _not_ghost, _ek_structural)) { #define COMPUTE_STRESS_ON_QUAD(type) computeStressOnQuad(); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(COMPUTE_STRESS_ON_QUAD); #undef COMPUTE_STRESS_ON_QUAD } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ bool StructuralMechanicsModel::allocateLumpedMassArray() { if (this->mass != nullptr) // Already allocated, so nothing to do. { return true; }; // now allocate it this->allocNodalField(this->mass, this->nb_degree_of_freedom, "lumped_mass"); return true; }; /* -------------------------------------------------------------------------- */ std::shared_ptr StructuralMechanicsModel::createNodalFieldBool( const std::string & field_name, const std::string & group_name, __attribute__((unused)) bool padding_flag) { std::map *> uint_nodal_fields; uint_nodal_fields["blocked_dofs"] = blocked_dofs.get(); return mesh.createNodalField(uint_nodal_fields[field_name], group_name); } /* -------------------------------------------------------------------------- */ std::shared_ptr StructuralMechanicsModel::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { UInt n; if (spatial_dimension == 2) { n = 2; } else { n = 3; } UInt padding_size = 0; if (padding_flag) { padding_size = 3; } if (field_name == "displacement") { return mesh.createStridedNodalField(displacement_rotation.get(), group_name, n, 0, padding_size); } if (field_name == "velocity") { return mesh.createStridedNodalField(velocity.get(), group_name, n, 0, padding_size); } if (field_name == "acceleration") { return mesh.createStridedNodalField(acceleration.get(), group_name, n, 0, padding_size); } if (field_name == "rotation") { return mesh.createStridedNodalField(displacement_rotation.get(), group_name, nb_degree_of_freedom - n, n, padding_size); } if (field_name == "force") { return mesh.createStridedNodalField(external_force.get(), group_name, n, 0, padding_size); } if (field_name == "external_force") { return mesh.createStridedNodalField(external_force.get(), group_name, n, 0, padding_size); } if (field_name == "momentum") { return mesh.createStridedNodalField(external_force.get(), group_name, nb_degree_of_freedom - n, n, padding_size); } if (field_name == "internal_force") { return mesh.createStridedNodalField(internal_force.get(), group_name, n, 0, padding_size); } if (field_name == "internal_momentum") { return mesh.createStridedNodalField(internal_force.get(), group_name, nb_degree_of_freedom - n, n, padding_size); } if (field_name == "mass") { AKANTU_DEBUG_ASSERT(this->mass.get() != nullptr, "The lumped mass matrix was not allocated."); return mesh.createStridedNodalField(this->mass.get(), group_name, n, 0, padding_size); } return nullptr; } /* -------------------------------------------------------------------------- */ std::shared_ptr StructuralMechanicsModel::createElementalField( const std::string & field_name, const std::string & group_name, bool /*unused*/, UInt spatial_dimension, ElementKind kind) { std::shared_ptr field; if (field_name == "element_index_by_material") { field = mesh.createElementalField( field_name, group_name, spatial_dimension, kind); } if (field_name == "stress") { ElementTypeMap nb_data_per_elem = this->mesh.getNbDataPerElem(stress); field = mesh.createElementalField( stress, group_name, this->spatial_dimension, kind, nb_data_per_elem); } return field; } /* -------------------------------------------------------------------------- */ /* Virtual methods from SolverCallback */ /* -------------------------------------------------------------------------- */ /// get the type of matrix needed MatrixType StructuralMechanicsModel::getMatrixType(const ID & /*id*/) const { return _symmetric; } /// callback to assemble a Matrix void StructuralMechanicsModel::assembleMatrix(const ID & id) { if (id == "K") { assembleStiffnessMatrix(); } else if (id == "M") { assembleMassMatrix(); } } /// callback to assemble a lumped Matrix void StructuralMechanicsModel::assembleLumpedMatrix(const ID & id) { if ("M" == id) { this->assembleLumpedMassMatrix(); } return; } /// callback to assemble the residual StructuralMechanicsModel::(rhs) void StructuralMechanicsModel::assembleResidual() { AKANTU_DEBUG_IN(); auto & dof_manager = getDOFManager(); assembleInternalForce(); // Ensures that the matrix are assembled. if (dof_manager.hasMatrix("K")) { this->assembleMatrix("K"); } if (dof_manager.hasMatrix("M")) { this->assembleMatrix("M"); } if (dof_manager.hasLumpedMatrix("M")) { this->assembleLumpedMassMatrix(); } /* This is essentially a summing up of forces * first the external forces are counted for and then stored inside the * residual. */ dof_manager.assembleToResidual("displacement", *external_force, 1); dof_manager.assembleToResidual("displacement", *internal_force, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleResidual(const ID & residual_part) { AKANTU_DEBUG_IN(); - // Ensures that the matrix are assembled. - if (dof_manager->hasMatrix("K")) { - this->assembleMatrix("K"); - } - if (dof_manager->hasMatrix("M")) { - this->assembleMatrix("M"); - } - if (dof_manager->hasLumpedMatrix("M")) { - this->assembleLumpedMassMatrix(); - } + auto & dof_manager = this->getDOFManager(); if ("external" == residual_part) { - this->getDOFManager().assembleToResidual("displacement", + dof_manager.assembleToResidual("displacement", *this->external_force, 1); AKANTU_DEBUG_OUT(); return; } if ("internal" == residual_part) { this->assembleInternalForce(); - this->getDOFManager().assembleToResidual("displacement", + dof_manager.assembleToResidual("displacement", *this->internal_force, 1); AKANTU_DEBUG_OUT(); return; } AKANTU_CUSTOM_EXCEPTION( debug::SolverCallbackResidualPartUnknown(residual_part)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Virtual methods from Model */ /* -------------------------------------------------------------------------- */ /// get some default values for derived classes std::tuple StructuralMechanicsModel::getDefaultSolverID(const AnalysisMethod & method) { switch (method) { case _static: { return std::make_tuple("static", TimeStepSolverType::_static); } case _implicit_dynamic: { return std::make_tuple("implicit", TimeStepSolverType::_dynamic); } case _explicit_lumped_mass: { // Taken from the solid mechanics part return std::make_tuple("explicit_lumped", TimeStepSolverType::_dynamic_lumped); } case _explicit_consistent_mass: { // Taken from the solid mechanics part return std::make_tuple("explicit", TimeStepSolverType::_dynamic); } default: std::cout << "UNKOWN." << std::endl; return std::make_tuple("unknown", TimeStepSolverType::_not_defined); } } /* ------------------------------------------------------------------------ */ ModelSolverOptions StructuralMechanicsModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_dynamic_lumped: { // Taken from the solid mechanic // part 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; // _linear; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } #if 1 case TimeStepSolverType::_dynamic: { // Copied from solid 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; } #else case TimeStepSolverType::_dynamic: { // Original options.non_linear_solver_type = NonLinearSolverType::_newton_raphson; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_trapezoidal_rule_2; options.solution_type["displacement"] = IntegrationScheme::_displacement; break; } #endif default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); } return options; } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleInternalForce() { internal_force->zero(); computeStresses(); for (auto type : mesh.elementTypes(_spatial_dimension = _all_dimensions, _element_kind = _ek_structural)) { assembleInternalForce(type, _not_ghost); // assembleInternalForce(type, _ghost); } } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleInternalForce(ElementType type, GhostType gt) { auto & fem = getFEEngine(); auto & sigma = stress(type, gt); auto ndof = getNbDegreeOfFreedom(type); auto nb_nodes = mesh.getNbNodesPerElement(type); auto ndof_per_elem = ndof * nb_nodes; Array BtSigma(fem.getNbIntegrationPoints(type) * mesh.getNbElement(type), ndof_per_elem, "BtSigma"); fem.computeBtD(sigma, BtSigma, type, gt); Array intBtSigma(0, ndof_per_elem, "intBtSigma"); fem.integrate(BtSigma, intBtSigma, ndof_per_elem, type, gt); getDOFManager().assembleElementalArrayLocalArray(intBtSigma, *internal_force, type, gt, -1.); } /* -------------------------------------------------------------------------- */ Real StructuralMechanicsModel::getKineticEnergy() { const UInt nb_nodes = mesh.getNbNodes(); const UInt nb_degree_of_freedom = this->nb_degree_of_freedom; Real ekin = 0.; // used to sum up energy (is divided by two at the very end) + //if mass matrix was not assembled, assemble it now + this->assembleMassMatrix(); + if (this->getDOFManager().hasLumpedMatrix("M")) { /* This code computes the kinetic energy for the case when the mass is * lumped. It is based on the solid mechanic equivalent. */ AKANTU_DEBUG_ASSERT(this->mass != nullptr, "The lumped mass is not allocated."); if (this->need_to_reassemble_lumpedMass) { this->assembleLumpedMatrix("M"); } /* Iterating over all nodes. * Important the velocity and mass also contains the rotational parts. * However, they can be handled in an uniform way. */ for (auto && data : zip(arange(nb_nodes), make_view(*this->velocity, nb_degree_of_freedom), make_view(*this->mass, nb_degree_of_freedom))) { const UInt n = std::get<0>(data); // This is the ID of the current node if (not mesh.isLocalOrMasterNode( n)) // Only handle the node if it belongs to us. { continue; } const auto & v = std::get<1>(data); // Get the velocity and mass of that node. const auto & m = std::get<2>(data); Real mv2 = 0.; // Contribution of this node. for (UInt i = 0; i < nb_degree_of_freedom; ++i) { /* In the solid mechanics part, only masses that are above a certain * value are considered. * However, the structural part, does not do this. */ const Real v_ = v(i); const Real m_ = m(i); mv2 += v_ * v_ * m_; } // end for(i): going through the components ekin += mv2; // add continution } // end for(n): iterating through all nodes } else if (this->getDOFManager().hasMatrix("M")) { /* Handle the case where no lumped mass is there. * This is basically the original code. */ if (this->need_to_reassemble_mass) { this->assembleMassMatrix(); } Array Mv(nb_nodes, nb_degree_of_freedom); this->getDOFManager().assembleMatMulVectToArray("displacement", "M", *this->velocity, Mv); for (auto && data : zip(arange(nb_nodes), make_view(Mv, nb_degree_of_freedom), make_view(*this->velocity, nb_degree_of_freedom))) { if (mesh.isLocalOrMasterNode(std::get<0>( data))) // only consider the node if we are blonging to it { ekin += std::get<2>(data).dot(std::get<1>(data)); } } } else { /* This is the case where no mass is present, for whatever reason, such as * the static case. We handle it specially be returning directly zero. * However, by doing that there will not be a syncronizing event as in the * other cases. Which is faster, but could be a problem in case the user * expects this. * * Another not is, that the solid mechanics part, would generate an error in * this clause. But, since the original implementation of the structural * part, did not do that, I, Philip, decided to refrain from that. However, * it is an option that should be considered. */ return 0.; } // Sum up across the comunicator mesh.getCommunicator().allReduce(ekin, SynchronizerOperation::_sum); return ekin / 2.; // finally divide the energy by two } /* -------------------------------------------------------------------------- */ Real StructuralMechanicsModel::getPotentialEnergy() { Real epot = 0.; UInt nb_nodes = mesh.getNbNodes(); - // if stiffness matrix is not assembled, do it - // as an alternative, gernate an error. - if (this->need_to_reassemble_stiffness) { - this->assembleStiffnessMatrix(); - }; + //if stiffness matrix is not assembled, do it + this->assembleStiffnessMatrix(); Array Ku(nb_nodes, nb_degree_of_freedom); this->getDOFManager().assembleMatMulVectToArray( "displacement", "K", *this->displacement_rotation, Ku); for (auto && data : zip(arange(nb_nodes), make_view(Ku, nb_degree_of_freedom), make_view(*this->displacement_rotation, nb_degree_of_freedom))) { epot += std::get<2>(data).dot(std::get<1>(data)) * static_cast(mesh.isLocalOrMasterNode(std::get<0>(data))); } mesh.getCommunicator().allReduce(epot, SynchronizerOperation::_sum); return epot / 2.; } /* -------------------------------------------------------------------------- */ Real StructuralMechanicsModel::getEnergy(const ID & energy) { if (energy == "kinetic") { return getKineticEnergy(); } if (energy == "potential") { return getPotentialEnergy(); } return 0; } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::computeForcesByLocalTractionArray( const Array & tractions, ElementType type) { AKANTU_DEBUG_IN(); auto nb_element = getFEEngine().getMesh().getNbElement(type); auto nb_nodes_per_element = getFEEngine().getMesh().getNbNodesPerElement(type); auto nb_quad = getFEEngine().getNbIntegrationPoints(type); // check dimension match AKANTU_DEBUG_ASSERT( Mesh::getSpatialDimension(type) == getFEEngine().getElementDimension(), "element type dimension does not match the dimension of boundaries : " << getFEEngine().getElementDimension() << " != " << Mesh::getSpatialDimension(type)); // check size of the vector AKANTU_DEBUG_ASSERT( tractions.size() == nb_quad * nb_element, "the size of the vector should be the total number of quadrature points"); // check number of components AKANTU_DEBUG_ASSERT(tractions.getNbComponent() == nb_degree_of_freedom, "the number of components should be the spatial " "dimension of the problem"); Array Ntbs(nb_element * nb_quad, nb_degree_of_freedom * nb_nodes_per_element); auto & fem = getFEEngine(); fem.computeNtb(tractions, Ntbs, type); // allocate the vector that will contain the integrated values auto name = id + std::to_string(type) + ":integral_boundary"; Array int_funct(nb_element, nb_degree_of_freedom * nb_nodes_per_element, name); // do the integration getFEEngine().integrate(Ntbs, int_funct, nb_degree_of_freedom * nb_nodes_per_element, type); // assemble the result into force vector getDOFManager().assembleElementalArrayLocalArray(int_funct, *external_force, type, _not_ghost, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::computeForcesByGlobalTractionArray( const Array & traction_global, ElementType type) { AKANTU_DEBUG_IN(); UInt nb_element = mesh.getNbElement(type); UInt nb_quad = getFEEngine().getNbIntegrationPoints(type); Array traction_local(nb_element * nb_quad, nb_degree_of_freedom, id + ":structuralmechanics:imposed_linear_load"); auto R_it = getFEEngineClass() .getShapeFunctions() .getRotations(type) .begin(nb_degree_of_freedom, nb_degree_of_freedom); auto Te_it = traction_global.begin(nb_degree_of_freedom); auto te_it = traction_local.begin(nb_degree_of_freedom); for (UInt e = 0; e < nb_element; ++e, ++R_it) { for (UInt q = 0; q < nb_quad; ++q, ++Te_it, ++te_it) { // turn the traction in the local referential te_it->template mul(*R_it, *Te_it); } } computeForcesByLocalTractionArray(traction_local, type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::afterSolveStep(bool converged) { if (converged) { assembleInternalForce(); } } } // namespace akantu