diff --git a/src/model/structural_mechanics/structural_mechanics_model.cc b/src/model/structural_mechanics/structural_mechanics_model.cc index 3098b1b22..de1ef7b7a 100644 --- a/src/model/structural_mechanics/structural_mechanics_model.cc +++ b/src/model/structural_mechanics/structural_mechanics_model.cc @@ -1,483 +1,482 @@ /** * @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: Wed Feb 21 2018 * * @brief Model implementation for Structural Mechanics elements * * * 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 "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" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER #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" #endif /* -------------------------------------------------------------------------- */ #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, const MemoryID & memory_id) : Model(mesh, ModelType::_structural_mechanics_model, dim, id, memory_id), f_m2a(1.0), stress("stress", id, memory_id), element_material("element_material", id, memory_id), set_ID("beam sets", id, memory_id), rotation_matrix("rotation_matices", id, memory_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(); } #ifdef AKANTU_USE_IOHELPER this->mesh.registerDumper("structural_mechanics_model", id, true); #endif 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); #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper().setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ /* 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); } } 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(); } /* -------------------------------------------------------------------------- */ 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; 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, group_name, n, 0, padding_size); } if (field_name == "velocity") { return mesh.createStridedNodalField(velocity, group_name, n, 0, padding_size); } if (field_name == "acceleration") { return mesh.createStridedNodalField(acceleration, group_name, n, 0, padding_size); } if (field_name == "rotation") { return mesh.createStridedNodalField(displacement_rotation, group_name, nb_degree_of_freedom - n, n, padding_size); } if (field_name == "force") { return mesh.createStridedNodalField(external_force, group_name, n, 0, padding_size); } if (field_name == "external_force") { return mesh.createStridedNodalField(external_force, group_name, n, 0, padding_size); } if (field_name == "momentum") { return mesh.createStridedNodalField( external_force, group_name, nb_degree_of_freedom - n, n, padding_size); } if (field_name == "internal_force") { return mesh.createStridedNodalField(internal_force, group_name, n, 0, padding_size); } if (field_name == "internal_momentum") { return mesh.createStridedNodalField( internal_force, group_name, nb_degree_of_freedom - n, n, 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*/) { 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*/) {} /// callback to assemble the residual StructuralMechanicsModel::(rhs) void StructuralMechanicsModel::assembleResidual() { AKANTU_DEBUG_IN(); auto & dof_manager = getDOFManager(); assembleInternalForce(); 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(); if ("external" == residual_part) { this->getDOFManager().assembleToResidual("displacement", *this->external_force, 1); AKANTU_DEBUG_OUT(); return; } if ("internal" == residual_part) { this->assembleInternalForce(); this->getDOFManager().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); } default: return std::make_tuple("unknown", TimeStepSolverType::_not_defined); } } /* ------------------------------------------------------------------------ */ ModelSolverOptions StructuralMechanicsModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case TimeStepSolverType::_static: { options.non_linear_solver_type = NonLinearSolverType::_linear; options.integration_scheme_type["displacement"] = IntegrationSchemeType::_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } case TimeStepSolverType::_dynamic: { 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; } /* -------------------------------------------------------------------------- */ 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); - BtSigma.resize(0); getDOFManager().assembleElementalArrayLocalArray(intBtSigma, *internal_force, type, gt, -1.); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/structural_mechanics/structural_mechanics_model.hh b/src/model/structural_mechanics/structural_mechanics_model.hh index f8a6db0db..48f0bbc68 100644 --- a/src/model/structural_mechanics/structural_mechanics_model.hh +++ b/src/model/structural_mechanics/structural_mechanics_model.hh @@ -1,351 +1,338 @@ /** * @file structural_mechanics_model.hh * * @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: Tue Feb 20 2018 * * @brief Particular implementation of the structural elements in the * StructuralMechanicsModel * * * 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 "aka_named_argument.hh" #include "boundary_condition.hh" #include "model.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_STRUCTURAL_MECHANICS_MODEL_HH_ #define AKANTU_STRUCTURAL_MECHANICS_MODEL_HH_ /* -------------------------------------------------------------------------- */ namespace akantu { class Material; class MaterialSelector; class DumperIOHelper; class NonLocalManager; template class IntegratorGauss; template class ShapeStructural; } // namespace akantu namespace akantu { struct StructuralMaterial { Real E{0}; Real A{1}; Real I{0}; Real Iz{0}; Real Iy{0}; Real GJ{0}; Real rho{0}; Real t{0}; Real nu{0}; }; class StructuralMechanicsModel : public Model { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: using MyFEEngineType = FEEngineTemplate; StructuralMechanicsModel(Mesh & mesh, UInt dim = _all_dimensions, const ID & id = "structural_mechanics_model", const MemoryID & memory_id = 0); ~StructuralMechanicsModel() override; /// Init full model void initFullImpl(const ModelOptions & options) override; /// Init boundary FEEngine void initFEEngineBoundary() override; /* ------------------------------------------------------------------------ */ /* Virtual methods from SolverCallback */ /* ------------------------------------------------------------------------ */ /// get the type of matrix needed MatrixType getMatrixType(const ID & matrix_id) override; /// callback to assemble a Matrix void assembleMatrix(const ID & matrix_id) override; /// callback to assemble a lumped Matrix void assembleLumpedMatrix(const ID & matrix_id) override; /// callback to assemble the residual (rhs) void assembleResidual() override; void assembleResidual(const ID & residual_part) override; bool canSplitResidual() override { return false; } /* ------------------------------------------------------------------------ */ /* Virtual methods from Model */ /* ------------------------------------------------------------------------ */ protected: /// get some default values for derived classes std::tuple getDefaultSolverID(const AnalysisMethod & method) override; ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; static UInt getNbDegreeOfFreedom(ElementType type); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ void initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) override; /// initialize the model void initModel() override; /// compute the stresses per elements void computeStresses(); /// compute the nodal forces void assembleInternalForce(); /// compute the nodal forces for an element type void assembleInternalForce(ElementType type, GhostType gt); /// assemble the stiffness matrix void assembleStiffnessMatrix(); /// assemble the mass matrix for consistent mass resolutions void assembleMassMatrix(); protected: /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMassMatrix(GhostType ghost_type); /// computes rho void computeRho(Array & rho, ElementType type, GhostType ghost_type); /// finish the computation of residual to solve in increment void updateResidualInternal(); /* ------------------------------------------------------------------------ */ private: template void assembleStiffnessMatrix(); template void computeStressOnQuad(); template void computeTangentModuli(Array & tangent_moduli); /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: std::shared_ptr createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; std::shared_ptr createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, UInt spatial_dimension, ElementKind kind) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// set the value of the time step void setTimeStep(Real time_step, const ID & solver_id = "") override; /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the StructuralMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement_rotation, Array &); /// get the StructuralMechanicsModel::velocity vector AKANTU_GET_MACRO(Velocity, *velocity, Array &); /// get the StructuralMechanicsModel::acceleration vector, updated /// by /// StructuralMechanicsModel::updateAcceleration AKANTU_GET_MACRO(Acceleration, *acceleration, Array &); /// get the StructuralMechanicsModel::external_force vector AKANTU_GET_MACRO(ExternalForce, *external_force, Array &); /// get the StructuralMechanicsModel::internal_force vector (boundary forces) AKANTU_GET_MACRO(InternalForce, *internal_force, Array &); /// get the StructuralMechanicsModel::boundary vector AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(RotationMatrix, rotation_matrix, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementMaterial, element_material, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Set_ID, set_ID, UInt); /** - * \brief This function adds the `StructuralMaterial` material to the list of materials managed by *this. + * \brief This function adds the `StructuralMaterial` material to the list of + * materials managed by *this. * - * It is important that this function might invalidate all references to structural materials, - * that were previously optained by `getMaterial()`. + * It is important that this function might invalidate all references to + * structural materials, that were previously optained by `getMaterial()`. * - * \param material The new material. + * \param material The new material. * - * \return The ID of the material that was added. + * \return The ID of the material that was added. * - * \note The return type is is new. + * \note The return type is is new. */ - UInt - addMaterial( - StructuralMaterial & material) - { - const auto matID = materials.size(); //ID of the material - materials.push_back(material); //add the material, might cause reallocation. - AKANTU_DEBUG_ASSERT(matID <= (::std::size_t)::std::numeric_limits::max(), "Can not represent the material ID"); + UInt addMaterial(StructuralMaterial & material) { + const auto matID = materials.size(); // ID of the material + materials.push_back(material); // add the material, might cause + // reallocation. + AKANTU_DEBUG_ASSERT(matID <= + (::std::size_t)::std::numeric_limits::max(), + "Can not represent the material ID"); return UInt(matID); } const StructuralMaterial & getMaterial(const Element & element) const { return materials[element_material(element)]; } - /** * \brief Returns the ith material of *this. * * \param i The ith material */ - const StructuralMaterial& - getMaterialByID( - UInt i) - const - noexcept(false) - { + const StructuralMaterial & getMaterialByID(UInt i) const noexcept(false) { return materials.at(i); } - /** * \brief Returns the number of the different materials inside *this. */ - UInt - getNbMaterials() - const - { - return materials.size(); - } - - + UInt getNbMaterials() const { return materials.size(); } /* ------------------------------------------------------------------------ */ /* Boundaries (structural_mechanics_model_boundary.cc) */ /* ------------------------------------------------------------------------ */ public: /// Compute Linear load function set in global axis template void computeForcesByGlobalTractionArray(const Array & tractions); /// Compute Linear load function set in local axis template void computeForcesByLocalTractionArray(const Array & tractions); /// compute force vector from a function(x,y,momentum) that describe stresses // template // void computeForcesFromFunction(BoundaryFunction in_function, // BoundaryFunctionType function_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// time step Real time_step; /// conversion coefficient form force/mass to acceleration Real f_m2a; /// displacements array Array * displacement_rotation{nullptr}; /// velocities array Array * velocity{nullptr}; /// accelerations array Array * acceleration{nullptr}; /// forces array Array * internal_force{nullptr}; /// forces array Array * external_force{nullptr}; /// lumped mass array Array * mass{nullptr}; /// boundaries array Array * blocked_dofs{nullptr}; /// stress array ElementTypeMapArray stress; ElementTypeMapArray element_material; // Define sets of beams ElementTypeMapArray set_ID; /// number of degre of freedom UInt nb_degree_of_freedom; // Rotation matrix ElementTypeMapArray rotation_matrix; // /// analysis method check the list in akantu::AnalysisMethod // AnalysisMethod method; /// flag defining if the increment must be computed or not bool increment_flag; bool need_to_reassemble_mass{true}; bool need_to_reassemble_stiffness{true}; /* ------------------------------------------------------------------------ */ std::vector materials; }; } // namespace akantu #include "structural_mechanics_model_inline_impl.hh" #endif /* AKANTU_STRUCTURAL_MECHANICS_MODEL_HH_ */ diff --git a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_2.cc b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_2.cc index e48ffd6b1..65c9f4dd6 100644 --- a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_2.cc +++ b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_2.cc @@ -1,112 +1,112 @@ /** * @file test_structural_mechanics_model_bernoulli_beam_2.cc * * @author Fabian Barras * @author Lucas Frerot * * @date creation: Fri Jul 15 2011 * @date last modification: Fri Feb 09 2018 * * @brief Computation of the analytical exemple 1.1 in the TGC vol 6 * * * 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 "test_structural_mechanics_model_fixture.hh" /* -------------------------------------------------------------------------- */ #include using namespace akantu; /* -------------------------------------------------------------------------- */ class TestStructBernoulli2 : public TestStructuralFixture> { using parent = TestStructuralFixture>; public: void addMaterials() override { mat.E = 3e10; mat.I = 0.0025; mat.A = 0.01; this->model->addMaterial(mat); mat.E = 3e10; mat.I = 0.00128; mat.A = 0.01; this->model->addMaterial(mat); } void assignMaterials() override { auto & materials = this->model->getElementMaterial(parent::type); materials(0) = 0; materials(1) = 1; } - void setDirichlets() override { + void setDirichletBCs() override { auto boundary = this->model->getBlockedDOFs().begin(parent::ndof); // clang-format off *boundary = {true, true, true}; ++boundary; *boundary = {false, true, false}; ++boundary; *boundary = {false, true, false}; ++boundary; // clang-format on } - void setNeumanns() override { + void setNeumannBCs() override { Real M = 3600; // Nm Real q = 6000; // kN/m Real L = 10; // m auto & forces = this->model->getExternalForce(); forces(2, 2) = -M; // moment on last node #if 1 // as long as integration is not available forces(0, 1) = -q * L / 2; forces(0, 2) = -q * L * L / 12; forces(1, 1) = -q * L / 2; forces(1, 2) = q * L * L / 12; #else auto & group = mesh.createElementGroup("lin_force"); group.add({type, 0, _not_ghost}); Vector lin_force = {0, q, 0}; // a linear force is not actually a *boundary* condition // it is equivalent to a volume force model.applyBC(BC::Neumann::FromSameDim(lin_force), group); #endif forces(2, 0) = mat.E * mat.A / 18; } protected: StructuralMaterial mat; }; /* -------------------------------------------------------------------------- */ TEST_F(TestStructBernoulli2, TestDisplacements) { model->solveStep(); auto d1 = model->getDisplacement()(1, 2); auto d2 = model->getDisplacement()(2, 2); auto d3 = model->getDisplacement()(1, 0); Real tol = Math::getTolerance(); EXPECT_NEAR(d1, 5.6 / 4800, tol); // first rotation EXPECT_NEAR(d2, -3.7 / 4800, tol); // second rotation EXPECT_NEAR(d3, 10. / 18, tol); // axial deformation } diff --git a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_3.cc b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_3.cc index e7d33e19b..4b334e128 100644 --- a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_3.cc +++ b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_3.cc @@ -1,96 +1,99 @@ /** * @file test_structural_mechanics_model_bernoulli_beam_3.cc * * @author Lucas Frerot * * @date creation: Fri Jul 15 2011 * @date last modification: Fri Feb 09 2018 * * @brief Computation of the analytical exemple 1.1 in the TGC vol 6 * * * 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 "test_structural_mechanics_model_fixture.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ using namespace akantu; class TestStructBernoulli3Static : public TestStructuralFixture> { using parent = TestStructuralFixture>; public: void readMesh(std::string filename) override { parent::readMesh(filename); + } + + void setNormals() override { auto &normals = this->mesh->getData("extra_normal", parent::type); normals(0, _z) = 1; normals(1, _z) = 1; } void addMaterials() override { StructuralMaterial mat; mat.E = 1; mat.Iz = 1; mat.Iy = 1; mat.A = 1; mat.GJ = 1; this->model->addMaterial(mat); } - void setDirichlets() override { + void setDirichletBCs() override { // Boundary conditions (blocking all DOFs of nodes 2 & 3) auto boundary = ++this->model->getBlockedDOFs().begin(parent::ndof); // clang-format off *boundary = {true, true, true, true, true, true}; ++boundary; *boundary = {true, true, true, true, true, true}; ++boundary; // clang-format on } - void setNeumanns() override { + void setNeumannBCs() override { // Forces Real P = 1; // N auto & forces = this->model->getExternalForce(); forces.zero(); forces(0, 2) = -P; // vertical force on first node } void assignMaterials() override { model->getElementMaterial(parent::type).set(0); } }; /* -------------------------------------------------------------------------- */ TEST_F(TestStructBernoulli3Static, TestDisplacements) { model->solveStep(); auto vz = model->getDisplacement()(0, 2); auto thy = model->getDisplacement()(0, 4); auto thx = model->getDisplacement()(0, 3); auto thz = model->getDisplacement()(0, 5); Real tol = Math::getTolerance(); EXPECT_NEAR(vz, -5. / 48, tol); EXPECT_NEAR(thy, -std::sqrt(2) / 8, tol); EXPECT_NEAR(thz, 0, tol); EXPECT_NEAR(thx, 0, tol); } diff --git a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_dynamics.cc b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_dynamics.cc index b28939ce3..c18ecc36a 100644 --- a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_dynamics.cc +++ b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_dynamics.cc @@ -1,177 +1,367 @@ /** * @file test_structural_mechanics_model_bernoulli_beam_dynamics.cc * * @author Sébastien Hartmann * * @date creation: Mon Jul 07 2014 * @date last modification: Wed Feb 03 2016 * * @brief Test for _bernouilli_beam in dynamic * * * Copyright (©) 2014-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 "test_structural_mechanics_model_fixture.hh" /* -------------------------------------------------------------------------- */ +#include "dof_manager.hh" #include "mesh_accessor.hh" #include "non_linear_solver_newton_raphson.hh" #include "structural_mechanics_model.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ using namespace akantu; /* -------------------------------------------------------------------------- */ static Real analytical_solution(Real time, Real L, Real rho, Real E, __attribute__((unused)) Real A, Real I, Real F) { Real omega = M_PI * M_PI / L / L * sqrt(E * I / rho); Real sum = 0.; UInt i = 5; for (UInt n = 1; n <= i; n += 2) { sum += (1. - cos(n * n * omega * time)) / pow(n, 4); } return 2. * F * pow(L, 3) / pow(M_PI, 4) / E / I * sum; } -class TestStructBernoulli3Dynamic - : public TestStructuralFixture> { - using parent = TestStructuralFixture>; +template +class TestStructBernoulliDynamic : public TestStructuralFixture { + using parent = TestStructuralFixture; public: - Real L{20}; - const UInt nb_element{10}; - UInt nb_nodes; + const UInt nb_element{40}; + const Real L{2}; + const Real le{L / nb_element}; + + const UInt nb_nodes{nb_element + 1}; + const Real F{1e4}; + StructuralMaterial mat; - const Real F{1.}; void readMesh(std::string /*filename*/) override { - nb_nodes = nb_element + 1; - auto length = L / nb_element; MeshAccessor mesh_accessor(*this->mesh); auto & nodes = mesh_accessor.getNodes(); nodes.resize(nb_nodes); - this->mesh->addConnectivityType(_bernoulli_beam_3); - auto & connectivities = mesh_accessor.getConnectivity(parent::type); - connectivities.resize(nb_element); - - this->mesh->getElementalData("extra_normal") - .initialize(*this->mesh, _element_kind = _ek_structural, - _nb_component = 3, _with_nb_element = true, - _default_value = 0.); - - auto & normals = this->mesh->getData("extra_normal", parent::type); - normals.resize(nb_element); + nodes.set(0.); - for (auto && data : enumerate(make_view(nodes, 3))) { + for (auto && data : enumerate(make_view(nodes, this->spatial_dimension))) { auto & node = std::get<1>(data); UInt i = std::get<0>(data); - node = {i * length, 0., 0.}; + node[_x] = i * le; } - for (auto && data : - enumerate(make_view(connectivities, 2), make_view(normals, 3))) { + this->mesh->addConnectivityType(parent::type); + auto & connectivities = mesh_accessor.getConnectivity(parent::type); + connectivities.resize(nb_element); + for (auto && data : enumerate(make_view(connectivities, 2))) { UInt i = std::get<0>(data); auto & connectivity = std::get<1>(data); - auto & normal = std::get<2>(data); connectivity = {i, i + 1}; - normal = {0., 0., 1.}; } mesh_accessor.makeReady(); } + void setNormals() override { + if (this->spatial_dimension != 3) { + return; + } + auto & normals = + this->mesh->template getData("extra_normal", parent::type); + normals.resize(nb_element); + + for (auto && normal : make_view(normals, this->spatial_dimension)) { + normal = {0., 0., 1.}; + } + } + AnalysisMethod getAnalysisMethod() const override { return _implicit_dynamic; } void addMaterials() override { this->mat.E = 1e9; - this->mat.rho = 1; + this->mat.rho = 10; + this->mat.I = 1; this->mat.Iz = 1; this->mat.Iy = 1; this->mat.A = 1; this->mat.GJ = 1; this->model->addMaterial(mat); } - void setDirichlets() override { - auto boundary = this->model->getBlockedDOFs().begin(parent::ndof); - // clang-format off - Vector boundary_b = boundary[0]; - Vector boundary_e = boundary[nb_nodes]; - boundary_b = {true, true, true, false, false, false}; - boundary_e = {false, true, true, false, false, false}; - // clang-format on + void setDirichletBCs() override { + auto & boundary = this->model->getBlockedDOFs(); + boundary.set(false); + + boundary(0, _x) = true; + boundary(0, _y) = true; + + boundary(nb_nodes - 1, _y) = true; + + if (this->spatial_dimension == 3) { + boundary(0, _z) = true; + boundary(nb_nodes - 1, _z) = true; + } } - void setNeumanns() override { - auto node_to_print = nb_nodes / 2 + 1; + void setNeumannBCs() override { + auto node_to_print = nb_nodes / 2; // Forces auto & forces = this->model->getExternalForce(); - forces(node_to_print - 1, _y) = F; + forces.zero(); + forces(node_to_print, _y) = F; } void assignMaterials() override { - model->getElementMaterial(parent::type).set(0); + this->model->getElementMaterial(parent::type).set(0); } }; -TEST_F(TestStructBernoulli3Dynamic, TestBeamOscilation) { +using beam_types = gtest_list_t, + element_type_t<_bernoulli_beam_3>>>; + +TYPED_TEST_SUITE(TestStructBernoulliDynamic, beam_types); + +template +void getElementMassMatrix(const StructuralMaterial & /*material*/, Real /*l*/, + Matrix & /*M*/) {} + +template +void getElementStifnessMatrix(const StructuralMaterial & /*material*/, + Real /*l*/, Matrix & /*M*/) {} + +template <> +void getElementMassMatrix>( + const StructuralMaterial & material, Real l, Matrix & M) { + auto A = material.A; + auto rho = material.rho; + // clang-format off + M = rho * A * l / 420. * Matrix({ + {140, 0, 0, 70, 0, 0}, + { 0, 156, 22*l, 0, 54, -13*l}, + { 0, 22*l, 4*l*l, 0, 13*l, -3*l*l}, + { 70, 0, 0, 140, 0, 0}, + { 0, 54, 13*l, 0, 156, -22*l}, + { 0,-13*l, -3*l*l, 0, -22*l, 4*l*l}}); + // clang-format on +} + +template <> +void getElementStifnessMatrix>( + const StructuralMaterial & material, Real l, Matrix & K) { + auto E = material.E; + auto A = material.A; + auto I = material.I; + + auto l_2 = l * l; + auto l_3 = l * l * l; + // clang-format off + K = Matrix({ + { E*A/l, 0, 0, -E*A/l, 0, 0}, + { 0, 12*E*I/l_3, 6*E*I/l_2, 0, -12*E*I/l_3, 6*E*I/l_2}, + { 0, 6*E*I/l_2, 4*E*I/l, 0, -6*E*I/l_2, 2*E*I/l}, + {-E*A/l, 0, 0, E*A/l, 0, 0}, + { 0, -12*E*I/l_3, -6*E*I/l_2, 0, 12*E*I/l_3, -6*E*I/l_2}, + { 0, 6*E*I/l_2, 2*E*I/l, 0, -6*E*I/l_2, 4*E*I/l}}); + // clang-format on +} + +template <> +void getElementMassMatrix>( + const StructuralMaterial & material, Real l, Matrix & M) { + auto A = material.A; + auto rho = material.rho; + // clang-format off + M = rho * A * l / 420. * Matrix({ + {140, 0, 0, 0, 0, 0, 70, 0, 0, 0, 0, 0}, + { 0, 156, 0, 0, 0, 22*l, 0, 54, 0, 0, 0, -13*l}, + { 0, 0, 156, 0, -22*l, 0, 0, 0, 54, 0, 13*l, 0}, + { 0, 0, 0, 140, 0, 0, 0, 0, 0, 70, 0, 0}, + { 0, 0, -22*l, 0, 4*l*l, 0, 0, 0, -13*l, 0, -3*l*l, 0}, + { 0, 22*l, 0, 0, 0, 4*l*l, 0, 13*l, 0, 0, 0, -3*l*l}, + { 70, 0, 0, 0, 0, 0, 140, 0, 0, 0, 0, 0}, + { 0, 54, 0, 0, 0, 13*l, 0, 156, 0, 0, 0, -22*l}, + { 0, 0, 54, 0, -13*l, 0, 0, 0, 156, 0, 22*l, 0}, + { 0, 0, 0, 70, 0, 0, 0, 0, 0, 140, 0, 0}, + { 0, 0, 13*l, 0, -3*l*l, 0, 0, 0, 22*l, 0, 4*l*l, 0}, + { 0, -13*l, 0, 0, 0, -3*l*l, 0, -22*l, 0, 0, 0, 4*l*l}}); + // clang-format on +} + +template <> +void getElementStifnessMatrix>( + const StructuralMaterial & material, Real l, Matrix & K) { + auto E = material.E; + auto A = material.A; + auto Iy = material.Iy; + auto Iz = material.Iz; + auto GJ = material.GJ; + + auto a1 = E * A / l; + auto b1 = 12 * E * Iz / l / l / l; + auto b2 = 6 * E * Iz / l / l; + auto b3 = 4 * E * Iz / l; + auto b4 = 2 * E * Iz / l; + + auto c1 = 12 * E * Iy / l / l / l; + auto c2 = 6 * E * Iy / l / l; + auto c3 = 4 * E * Iy / l; + auto c4 = 2 * E * Iy / l; + + auto d1 = GJ / l; + + // clang-format off + K = Matrix({ + { a1, 0, 0, 0, 0, 0, -a1, 0, 0, 0, 0, 0}, + { 0, b1, 0, 0, 0, b2, 0, -b1, 0, 0, 0, b2}, + { 0, 0, c1, 0, -c2, 0, 0, 0, -c1, 0, -c2, 0}, + { 0, 0, 0, d1, 0, 0, 0, 0, 0, -d1, 0, 0}, + { 0, 0, -c2, 0, c3, 0, 0, 0, c2, 0, c4, 0}, + { 0, b2, 0, 0, 0, b3, 0, -b2, 0, 0, 0, b4}, + { -a1, 0, 0, 0, 0, 0, a1, 0, 0, 0, 0, 0}, + { 0, -b1, 0, 0, 0, -b2, 0, b1, 0, 0, 0, -b2}, + { 0, 0, -c1, 0, c2, 0, 0, 0, c1, 0, c2, 0}, + { 0, 0, 0, -d1, 0, 0, 0, 0, 0, d1, 0, 0}, + { 0, 0, -c2, 0, c4, 0, 0, 0, c2, 0, c3, 0}, + { 0, b2, 0, 0, 0, b4, 0, -b2, 0, 0, 0, b3}}); + // clang-format on +} + +TYPED_TEST(TestStructBernoulliDynamic, TestBeamMatrices) { + this->model->assembleMatrix("M"); + this->model->assembleMatrix("K"); + + const auto & K = this->model->getDOFManager().getMatrix("K"); + const auto & M = this->model->getDOFManager().getMatrix("M"); + + Matrix Ka(this->nb_nodes * this->ndof, this->nb_nodes * this->ndof, 0.); + Matrix Ma(this->nb_nodes * this->ndof, this->nb_nodes * this->ndof, 0.); + + Matrix Ke(this->ndof * 2, this->ndof * 2); + Matrix Me(this->ndof * 2, this->ndof * 2); + + getElementMassMatrix(this->mat, this->le, Me); + getElementStifnessMatrix(this->mat, this->le, Ke); + + auto assemble = [&](auto && nodes, auto && M, auto && Me) { + auto n1 = nodes[0]; + auto n2 = nodes[1]; + for (auto i : arange(this->ndof)) { + for (auto j : arange(this->ndof)) { + M(n1 * this->ndof + i, n1 * this->ndof + j) += Me(i, j); + M(n2 * this->ndof + i, n2 * this->ndof + j) += + Me(this->ndof + i, this->ndof + j); + M(n1 * this->ndof + i, n2 * this->ndof + j) += Me(i, this->ndof + j); + M(n2 * this->ndof + i, n1 * this->ndof + j) += Me(this->ndof + i, j); + } + } + }; + + auto && connectivities = this->mesh->getConnectivity(this->type); + + for (auto && connectivity : make_view(connectivities, 2)) { + assemble(connectivity, Ka, Ke); + assemble(connectivity, Ma, Me); + } + + auto tol = 1e-13; + + auto Ka_max = Ka.template norm(); + auto Ma_max = Ma.template norm(); + + for (auto i : arange(Ka.rows())) { + for (auto j : arange(Ka.cols())) { + EXPECT_NEAR(Ka(i, j), K(i, j), tol * Ka_max); + EXPECT_NEAR(Ma(i, j), M(i, j), tol * Ma_max); + } + } +} + +TYPED_TEST(TestStructBernoulliDynamic, TestBeamOscilation) { + Real time_step = 1e-6; + this->model->setTimeStep(time_step); + auto & solver = this->model->getNonLinearSolver(); solver.set("max_iterations", 100); solver.set("threshold", 1e-8); solver.set("convergence_type", SolveConvergenceCriteria::_solution); - auto node_to_print = this->nb_nodes / 2 + 1; - auto & d = this->model->getDisplacement()(node_to_print, 1); - - Real time_step = 1e-5; - this->model->setTimeStep(time_step); - + auto node_to_print = this->nb_nodes / 2; + auto & d = this->model->getDisplacement()(node_to_print, _y); std::ofstream pos; - pos.open("position.csv"); + std::string filename = "position" + std::to_string(this->type) + ".csv"; + pos.open(filename); if (not pos.good()) { AKANTU_ERROR("Cannot open file \"position.csv\""); } pos << "id,time,position,solution" << std::endl; - Real tol = 1e-10; +//#define debug +#ifdef debug + this->model->addDumpFieldVector("displacement"); + this->model->addDumpField("blocked_dofs"); + this->model->addDumpFieldVector("external_force"); + this->model->addDumpFieldVector("internal_force"); + this->model->addDumpFieldVector("acceleration"); + this->model->addDumpFieldVector("velocity"); + this->model->dump(); +#endif + + this->model->getDisplacement().set(0.); + + Real tol = 1e-6; Real time = 0.; - for (UInt s = 1; time < 0.1606; ++s) { + for (UInt s = 1; s < 300; ++s) { EXPECT_NO_THROW(this->model->solveStep()); time = s * time_step; auto da = analytical_solution(time, this->L, this->mat.rho, this->mat.E, this->mat.A, this->mat.Iy, this->F); - pos << s << "," << time << "," << d << "," << da <model->dump(); +#endif + EXPECT_NEAR(d, da, tol); } + } diff --git a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_discrete_kirchhoff_triangle_18.cc b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_discrete_kirchhoff_triangle_18.cc index d254434cf..adbf05c80 100644 --- a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_discrete_kirchhoff_triangle_18.cc +++ b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_discrete_kirchhoff_triangle_18.cc @@ -1,111 +1,111 @@ /** * @file test_structural_mechanics_model_discrete_kirchhoff_triangle_18.cc * * @author Fabian Barras * @author Lucas Frerot * * @date creation: Fri Jul 15 2011 * @date last modification: Wed Feb 21 2018 * * @brief Computation of the analytical exemple 1.1 in the TGC vol 6 * * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "sparse_matrix.hh" #include "test_structural_mechanics_model_fixture.hh" /* -------------------------------------------------------------------------- */ #include using namespace akantu; /* -------------------------------------------------------------------------- */ class TestStructDKT18 : public TestStructuralFixture< element_type_t<_discrete_kirchhoff_triangle_18>> { using parent = TestStructuralFixture>; public: void addMaterials() override { mat.E = 1; mat.t = 1; mat.nu = 0.3; this->model->addMaterial(mat); } void assignMaterials() override { auto & materials = this->model->getElementMaterial(parent::type); std::fill(materials.begin(), materials.end(), 0); } - void setDirichlets() override { + void setDirichletBCs() override { this->model->getBlockedDOFs().set(true); auto center_node = this->model->getBlockedDOFs().end(parent::ndof) - 1; *center_node = {false, false, false, false, false, true}; this->model->getDisplacement().zero(); auto disp = ++this->model->getDisplacement().begin(parent::ndof); // Displacement field from Batoz Vol. 2 p. 392 // with theta to beta correction from discrete Kirchhoff condition // see also the master thesis of Michael Lozano // clang-format off // This displacement field tests membrane and bending modes *disp = {40, 20, -800 , -20, 40, 0}; ++disp; *disp = {50, 40, -1400, -40, 50, 0}; ++disp; *disp = {10, 20, -200 , -20, 10, 0}; ++disp; // This displacement tests the bending mode // *disp = {0, 0, -800 , -20, 40, 0}; ++disp; // *disp = {0, 0, -1400, -40, 50, 0}; ++disp; // *disp = {0, 0, -200 , -20, 10, 0}; ++disp; // This displacement tests the membrane mode // *disp = {40, 20, 0, 0, 0, 0}; ++disp; // *disp = {50, 40, 0, 0, 0, 0}; ++disp; // *disp = {10, 20, 0, 0, 0, 0}; ++disp; // clang-format on } - void setNeumanns() override {} + void setNeumannBCs() override {} protected: StructuralMaterial mat; }; /* -------------------------------------------------------------------------- */ // Batoz Vol 2. patch test, ISBN 2-86601-259-3 TEST_F(TestStructDKT18, TestDisplacements) { model->solveStep(); Vector solution = {22.5, 22.5, -337.5, -22.5, 22.5, 0}; auto nb_nodes = this->model->getDisplacement().size(); Vector center_node_disp = model->getDisplacement().begin(solution.size())[nb_nodes - 1]; auto error = solution - center_node_disp; EXPECT_NEAR(error.norm(), 0., 1e-12); } diff --git a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_fixture.hh b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_fixture.hh index a56601608..14b48d567 100644 --- a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_fixture.hh +++ b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_fixture.hh @@ -1,115 +1,115 @@ /** * @file test_structural_mechanics_model_fixture.hh * * @author Lucas Frerot * * @date creation: Tue Nov 14 2017 * @date last modification: Fri Feb 09 2018 * * @brief Main test for structural model * * * 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 "element_class_structural.hh" #include "structural_mechanics_model.hh" #include "test_gtest_utils.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_TEST_STRUCTURAL_MECHANICS_MODEL_FIXTURE_HH_ #define AKANTU_TEST_STRUCTURAL_MECHANICS_MODEL_FIXTURE_HH_ using namespace akantu; template class TestStructuralFixture : public ::testing::Test { public: static constexpr const ElementType type = type_::value; static constexpr const size_t spatial_dimension = ElementClass::getSpatialDimension(); static const UInt ndof = ElementClass::getNbDegreeOfFreedom(); void SetUp() override { const auto spatial_dimension = this->spatial_dimension; mesh = std::make_unique(spatial_dimension); readMesh(makeMeshName()); std::stringstream element_type; element_type << this->type; model = std::make_unique(*mesh, _all_dimensions, element_type.str()); + setNormals(); initModel(); } virtual void initModel() { this->addMaterials(); auto method = getAnalysisMethod(); + this->model->initFull(_analysis_method = method); + this->assignMaterials(); - this->setDirichlets(); - this->setNeumanns(); + + this->setDirichletBCs(); + this->setNeumannBCs(); } virtual AnalysisMethod getAnalysisMethod() const { return _static; } virtual void readMesh(std::string filename) { mesh->read(filename, _miot_gmsh_struct); - mesh->getElementalData("extra_normal") - .initialize(*this->mesh, _element_kind = _ek_structural, - _nb_component = spatial_dimension, _with_nb_element = true, - _default_value = 0.); } virtual std::string makeMeshName() { std::stringstream element_type; element_type << type; SCOPED_TRACE(element_type.str().c_str()); return element_type.str() + ".msh"; } void TearDown() override { model.reset(nullptr); mesh.reset(nullptr); } virtual void addMaterials() = 0; virtual void assignMaterials() = 0; - virtual void setDirichlets() = 0; - virtual void setNeumanns() = 0; + virtual void setDirichletBCs() = 0; + virtual void setNeumannBCs() = 0; + virtual void setNormals() {} protected: std::unique_ptr mesh; std::unique_ptr model; }; template constexpr ElementType TestStructuralFixture::type; template constexpr size_t TestStructuralFixture::spatial_dimension; template const UInt TestStructuralFixture::ndof; -// using types = gtest_list_t; +using structural_types = gtest_list_t; -// TYPED_TEST_SUITE(TestStructuralFixture, types); #endif /* AKANTU_TEST_STRUCTURAL_MECHANICS_MODEL_FIXTURE_HH_ */