diff --git a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.hh b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.hh index a6c3088ff..8ff027c8f 100644 --- a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.hh +++ b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement.hh @@ -1,213 +1,210 @@ /** * @file material_reinforcement.hh * * @author Lucas Frerot * * @date creation: Fri Mar 13 2015 * @date last modification: Tue Nov 24 2015 * * @brief Reinforcement material * * @section LICENSE * * Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory * (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_REINFORCEMENT_HH__ #define __AKANTU_MATERIAL_REINFORCEMENT_HH__ #include "aka_common.hh" #include "material.hh" #include "embedded_interface_model.hh" #include "embedded_internal_field.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /** * @brief Material used to represent embedded reinforcements * * This class is used for computing the reinforcement stiffness matrix * along with the reinforcement residual. Room is made for constitutive law, * but actual use of contitutive laws is made in MaterialReinforcementTemplate. * * Be careful with the dimensions in this class : * - this->spatial_dimension is always 1 * - the template parameter dim is the dimension of the problem */ template class MaterialReinforcement : public Mat { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /// Constructor MaterialReinforcement(EmbeddedInterfaceModel & model, const ID & id = ""); /// Destructor ~MaterialReinforcement() override; protected: void initialize(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// Init the material void initMaterial() override; /// Init the filters for background elements void initFilters(); /// Init the background shape derivatives void initBackgroundShapeDerivatives(); /// Init the cosine matrices void initDirectingCosines(); /// Assemble stiffness matrix void assembleStiffnessMatrix(GhostType ghost_type) override; /// Compute all the stresses ! void computeAllStresses(GhostType ghost_type) override; /// Compute energy Real getEnergy(const std::string & id) override; /// Assemble the residual of one type of element (typically _segment_2) void assembleInternalForces(GhostType ghost_type) override; /* ------------------------------------------------------------------------ */ /* Protected methods */ /* ------------------------------------------------------------------------ */ protected: /// Allocate the background shape derivatives void allocBackgroundShapeDerivatives(); /// Compute the directing cosines matrix for one element type void computeDirectingCosines(const ElementType & type, GhostType ghost_type); /// Compute the directing cosines matrix on quadrature points. inline void computeDirectingCosinesOnQuad(const Matrix & nodes, Matrix & cosines); /// Add the prestress to the computed stress void addPrestress(const ElementType & type, GhostType ghost_type); /// Compute displacement gradient in reinforcement void computeGradU(const ElementType & interface_type, GhostType ghost_type); /// Assemble the stiffness matrix for an element type (typically _segment_2) void assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type); /// Assemble the stiffness matrix for background & interface types void assembleStiffnessMatrixInterface(const ElementType & interface_type, const ElementType & background_type, GhostType ghost_type); /// Compute the background shape derivatives for a type void computeBackgroundShapeDerivatives(const ElementType & type, GhostType ghost_type); /// Compute the background shape derivatives for a type pair void computeBackgroundShapeDerivatives(const ElementType & interface_type, const ElementType & bg_type, GhostType ghost_type, const Array & filter); /// Filter elements crossed by interface of a type void filterInterfaceBackgroundElements(Array & foreground, Array & background, const ElementType & type, const ElementType & interface_type, GhostType ghost_type); /// Assemble the residual of one type of element (typically _segment_2) void assembleInternalForces(const ElementType & type, GhostType ghost_type); /// Assemble the residual for a pair of elements void assembleInternalForcesInterface(const ElementType & interface_type, const ElementType & background_type, GhostType ghost_type); // TODO figure out why voigt size is 4 in 2D inline void stressTensorToVoigtVector(const Matrix & tensor, Vector & vector); inline void strainTensorToVoigtVector(const Matrix & tensor, Vector & vector); /// Get background filter Array & getBackgroundFilter(const ElementType & fg_type, const ElementType & bg_type, GhostType ghost_type) { return (*background_filter(fg_type, ghost_type))(bg_type, ghost_type); } /// Get foreground filter Array & getForegroundFilter(const ElementType & fg_type, const ElementType & bg_type, GhostType ghost_type) { return (*foreground_filter(fg_type, ghost_type))(bg_type, ghost_type); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// Embedded model EmbeddedInterfaceModel & emodel; - /// Stress in the reinforcement - InternalField stress_embedded; - /// Gradu of concrete on reinforcement InternalField gradu_embedded; /// C matrix on quad InternalField directing_cosines; /// Prestress on quad InternalField pre_stress; /// Cross-sectional area Real area; template using CrossMap = ElementTypeMap>>; /// Background mesh shape derivatives CrossMap shape_derivatives; /// Foreground mesh filter (contains segment ids) CrossMap foreground_filter; /// Background element filter (contains bg ids) CrossMap background_filter; }; } // akantu #include "material_reinforcement_tmpl.hh" #endif // __AKANTU_MATERIAL_REINFORCEMENT_HH__ diff --git a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_tmpl.hh b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_tmpl.hh index 17e60356b..09caa4037 100644 --- a/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_tmpl.hh +++ b/src/model/solid_mechanics/materials/material_embedded/material_reinforcement_tmpl.hh @@ -1,822 +1,818 @@ /** * @file material_reinforcement_tmpl.hh * * @author Lucas Frerot * * @date creation: Thu Feb 1 2018 * * @brief Reinforcement material * * @section LICENSE * * Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory * (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_voigthelper.hh" #include "material_reinforcement.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialReinforcement::MaterialReinforcement( EmbeddedInterfaceModel & model, const ID & id) : Mat(model, 1, model.getInterfaceMesh(), model.getFEEngine("EmbeddedInterfaceFEEngine"), id), emodel(model), - stress_embedded("stress_embedded", *this, 1, - model.getFEEngine("EmbeddedInterfaceFEEngine"), - this->element_filter), gradu_embedded("gradu_embedded", *this, 1, model.getFEEngine("EmbeddedInterfaceFEEngine"), this->element_filter), directing_cosines("directing_cosines", *this, 1, model.getFEEngine("EmbeddedInterfaceFEEngine"), this->element_filter), pre_stress("pre_stress", *this, 1, model.getFEEngine("EmbeddedInterfaceFEEngine"), this->element_filter), area(1.0), shape_derivatives() { AKANTU_DEBUG_IN(); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::initialize() { AKANTU_DEBUG_IN(); this->registerParam("area", area, _pat_parsable | _pat_modifiable, "Reinforcement cross-sectional area"); this->registerParam("pre_stress", pre_stress, _pat_parsable | _pat_modifiable, "Uniform pre-stress"); // this->unregisterInternal(this->stress); // Fool the AvgHomogenizingFunctor // stress.initialize(dim * dim); // Reallocate the element filter this->element_filter.initialize(this->emodel.getInterfaceMesh(), _spatial_dimension = 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template MaterialReinforcement::~MaterialReinforcement() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::initMaterial() { Mat::initMaterial(); - stress_embedded.initialize(dim * dim); gradu_embedded.initialize(dim * dim); pre_stress.initialize(1); /// We initialise the stuff that is not going to change during the simulation this->initFilters(); this->allocBackgroundShapeDerivatives(); this->initBackgroundShapeDerivatives(); this->initDirectingCosines(); } /* -------------------------------------------------------------------------- */ namespace detail { class FilterInitializer : public MeshElementTypeMapArrayInializer { public: FilterInitializer(EmbeddedInterfaceModel & emodel, const GhostType & ghost_type) : MeshElementTypeMapArrayInializer(emodel.getMesh(), 1, emodel.getSpatialDimension(), ghost_type, _ek_regular) {} UInt size(const ElementType & /*bgtype*/) const override { return 0; } }; } /* -------------------------------------------------------------------------- */ /// Initialize the filter for background elements template void MaterialReinforcement::initFilters() { for (auto gt : ghost_types) { for (auto && type : emodel.getInterfaceMesh().elementTypes(1, gt)) { std::string shaped_id = "filter"; if (gt == _ghost) shaped_id += ":ghost"; auto & background = background_filter(std::make_unique>( "bg_" + shaped_id, this->name), type, gt); auto & foreground = foreground_filter( std::make_unique>(shaped_id, this->name), type, gt); foreground->initialize(detail::FilterInitializer(emodel, gt), 0, true); background->initialize(detail::FilterInitializer(emodel, gt), 0, true); // Computing filters for (auto && bg_type : background->elementTypes(dim, gt)) { filterInterfaceBackgroundElements( (*foreground)(bg_type), (*background)(bg_type), bg_type, type, gt); } } } } /* -------------------------------------------------------------------------- */ /// Construct a filter for a (interface_type, background_type) pair template void MaterialReinforcement::filterInterfaceBackgroundElements( Array & foreground, Array & background, const ElementType & type, const ElementType & interface_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); foreground.resize(0); background.resize(0); Array & elements = emodel.getInterfaceAssociatedElements(interface_type, ghost_type); Array & elem_filter = this->element_filter(interface_type, ghost_type); for (auto & elem_id : elem_filter) { Element & elem = elements(elem_id); if (elem.type == type) { background.push_back(elem.element); foreground.push_back(elem_id); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ namespace detail { class BackgroundShapeDInitializer : public ElementTypeMapArrayInializer { public: BackgroundShapeDInitializer(UInt spatial_dimension, FEEngine & engine, const ElementType & foreground_type, ElementTypeMapArray & filter, const GhostType & ghost_type) : ElementTypeMapArrayInializer(spatial_dimension, 0, ghost_type, _ek_regular) { auto nb_quad = engine.getNbIntegrationPoints(foreground_type); // Counting how many background elements are affected by elements of // interface_type for (auto type : filter.elementTypes(this->spatial_dimension)) { // Inserting size array_size_per_bg_type(filter(type).size() * nb_quad, type, this->ghost_type); } } auto elementTypes() const -> decltype(auto) { return array_size_per_bg_type.elementTypes(); } UInt size(const ElementType & bgtype) const { return array_size_per_bg_type(bgtype, this->ghost_type); } UInt nbComponent(const ElementType & bgtype) const { return ShapeFunctions::getShapeDerivativesSize(bgtype); } protected: ElementTypeMap array_size_per_bg_type; }; } /** * Background shape derivatives need to be stored per background element * types but also per embedded element type, which is why they are stored * in an ElementTypeMap *>. The outer ElementTypeMap * refers to the embedded types, and the inner refers to the background types. */ template void MaterialReinforcement::allocBackgroundShapeDerivatives() { AKANTU_DEBUG_IN(); for (auto gt : ghost_types) { for (auto && type : emodel.getInterfaceMesh().elementTypes(1, gt)) { std::string shaped_id = "embedded_shape_derivatives"; if (gt == _ghost) shaped_id += ":ghost"; auto & shaped_etma = shape_derivatives( std::make_unique>(shaped_id, this->name), type, gt); shaped_etma->initialize( detail::BackgroundShapeDInitializer( emodel.getSpatialDimension(), emodel.getFEEngine("EmbeddedInterfaceFEEngine"), type, *background_filter(type, gt), gt), 0, true); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::initBackgroundShapeDerivatives() { AKANTU_DEBUG_IN(); for (auto interface_type : this->element_filter.elementTypes(this->spatial_dimension)) { for (auto type : background_filter(interface_type)->elementTypes(dim)) { computeBackgroundShapeDerivatives(interface_type, type, _not_ghost, this->element_filter(interface_type)); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::computeBackgroundShapeDerivatives( const ElementType & interface_type, const ElementType & bg_type, GhostType ghost_type, const Array & filter) { auto & interface_engine = emodel.getFEEngine("EmbeddedInterfaceFEEngine"); auto & engine = emodel.getFEEngine(); auto & interface_mesh = emodel.getInterfaceMesh(); const auto nb_nodes_elem_bg = Mesh::getNbNodesPerElement(bg_type); // const auto nb_strss = VoigtHelper::size; const auto nb_quads_per_elem = interface_engine.getNbIntegrationPoints(interface_type); Array quad_pos(0, dim, "interface_quad_pos"); interface_engine.interpolateOnIntegrationPoints(interface_mesh.getNodes(), quad_pos, dim, interface_type, ghost_type, filter); auto & background_shapesd = (*shape_derivatives(interface_type, ghost_type))(bg_type, ghost_type); auto & background_elements = (*background_filter(interface_type, ghost_type))(bg_type, ghost_type); auto & foreground_elements = (*foreground_filter(interface_type, ghost_type))(bg_type, ghost_type); auto shapesd_begin = background_shapesd.begin(dim, nb_nodes_elem_bg, nb_quads_per_elem); auto quad_begin = quad_pos.begin(dim, nb_quads_per_elem); for (auto && tuple : zip(background_elements, foreground_elements)) { UInt bg = std::get<0>(tuple), fg = std::get<1>(tuple); for (UInt i = 0; i < nb_quads_per_elem; ++i) { Matrix shapesd = Tensor3(shapesd_begin[fg])(i); Vector quads = Matrix(quad_begin[fg])(i); engine.computeShapeDerivatives(quads, bg, bg_type, shapesd, ghost_type); } } } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::initDirectingCosines() { AKANTU_DEBUG_IN(); Mesh & mesh = emodel.getInterfaceMesh(); Mesh::type_iterator type_it = mesh.firstType(1, _not_ghost); Mesh::type_iterator type_end = mesh.lastType(1, _not_ghost); const UInt voigt_size = VoigtHelper::size; directing_cosines.initialize(voigt_size); for (; type_it != type_end; ++type_it) { computeDirectingCosines(*type_it, _not_ghost); // computeDirectingCosines(*type_it, _ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::assembleStiffnessMatrix( GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh & interface_mesh = emodel.getInterfaceMesh(); Mesh::type_iterator type_it = interface_mesh.firstType(1, _not_ghost); Mesh::type_iterator type_end = interface_mesh.lastType(1, _not_ghost); for (; type_it != type_end; ++type_it) { assembleStiffnessMatrix(*type_it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::assembleInternalForces( GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh & interface_mesh = emodel.getInterfaceMesh(); Mesh::type_iterator type_it = interface_mesh.firstType(1, _not_ghost); Mesh::type_iterator type_end = interface_mesh.lastType(1, _not_ghost); for (; type_it != type_end; ++type_it) { this->assembleInternalForces(*type_it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::computeAllStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh::type_iterator it = emodel.getInterfaceMesh().firstType(); Mesh::type_iterator last_type = emodel.getInterfaceMesh().lastType(); for (; it != last_type; ++it) { computeGradU(*it, ghost_type); this->computeStress(*it, ghost_type); addPrestress(*it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::addPrestress(const ElementType & type, GhostType ghost_type) { auto & stress = this->stress(type, ghost_type); auto & sigma_p = this->pre_stress(type, ghost_type); for (auto && tuple : zip(stress, sigma_p)) { std::get<0>(tuple) += std::get<1>(tuple); } } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::assembleInternalForces( const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh & mesh = emodel.getMesh(); Mesh::type_iterator type_it = mesh.firstType(dim, ghost_type); Mesh::type_iterator type_end = mesh.lastType(dim, ghost_type); for (; type_it != type_end; ++type_it) { assembleInternalForcesInterface(type, *type_it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Computes and assemble the residual. Residual in reinforcement is computed as: * * \f[ * \vec{r} = A_s \int_S{\mathbf{B}^T\mathbf{C}^T \vec{\sigma_s}\,\mathrm{d}s} * \f] */ template void MaterialReinforcement::assembleInternalForcesInterface( const ElementType & interface_type, const ElementType & background_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt voigt_size = VoigtHelper::size; FEEngine & interface_engine = emodel.getFEEngine("EmbeddedInterfaceFEEngine"); Array & elem_filter = this->element_filter(interface_type, ghost_type); UInt nodes_per_background_e = Mesh::getNbNodesPerElement(background_type); UInt nb_quadrature_points = interface_engine.getNbIntegrationPoints(interface_type, ghost_type); UInt nb_element = elem_filter.size(); UInt back_dof = dim * nodes_per_background_e; Array & shapesd = (*shape_derivatives(interface_type, ghost_type))( background_type, ghost_type); Array integrant(nb_quadrature_points * nb_element, back_dof, "integrant"); Array::vector_iterator integrant_it = integrant.begin(back_dof); Array::vector_iterator integrant_end = integrant.end(back_dof); Array::matrix_iterator B_it = shapesd.begin(dim, nodes_per_background_e); Array::vector_iterator C_it = directing_cosines(interface_type, ghost_type).begin(voigt_size); auto sigma_it = this->stress(interface_type, ghost_type).begin(); Matrix Bvoigt(voigt_size, back_dof); for (; integrant_it != integrant_end; ++integrant_it, ++B_it, ++C_it, ++sigma_it) { VoigtHelper::transferBMatrixToSymVoigtBMatrix(*B_it, Bvoigt, nodes_per_background_e); Vector & C = *C_it; Vector & BtCt_sigma = *integrant_it; BtCt_sigma.mul(Bvoigt, C); BtCt_sigma *= *sigma_it * area; } Array residual_interface(nb_element, back_dof, "residual_interface"); interface_engine.integrate(integrant, residual_interface, back_dof, interface_type, ghost_type, elem_filter); integrant.resize(0); Array background_filter(nb_element, 1, "background_filter"); auto & filter = getBackgroundFilter(interface_type, background_type, ghost_type); emodel.getDOFManager().assembleElementalArrayLocalArray( residual_interface, emodel.getInternalForce(), background_type, ghost_type, -1., filter); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::computeDirectingCosines( const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh & interface_mesh = emodel.getInterfaceMesh(); const UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const UInt steel_dof = dim * nb_nodes_per_element; const UInt voigt_size = VoigtHelper::size; const UInt nb_quad_points = emodel.getFEEngine("EmbeddedInterfaceFEEngine") .getNbIntegrationPoints(type, ghost_type); Array node_coordinates(this->element_filter(type, ghost_type).size(), steel_dof); this->emodel.getFEEngine().template extractNodalToElementField( interface_mesh, interface_mesh.getNodes(), node_coordinates, type, ghost_type, this->element_filter(type, ghost_type)); Array::matrix_iterator directing_cosines_it = directing_cosines(type, ghost_type).begin(1, voigt_size); Array::matrix_iterator node_coordinates_it = node_coordinates.begin(dim, nb_nodes_per_element); Array::matrix_iterator node_coordinates_end = node_coordinates.end(dim, nb_nodes_per_element); for (; node_coordinates_it != node_coordinates_end; ++node_coordinates_it) { for (UInt i = 0; i < nb_quad_points; i++, ++directing_cosines_it) { Matrix & nodes = *node_coordinates_it; Matrix & cosines = *directing_cosines_it; computeDirectingCosinesOnQuad(nodes, cosines); } } // Mauro: the directing_cosines internal is defined on the quadrature points // of each element AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::assembleStiffnessMatrix( const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh & mesh = emodel.getMesh(); Mesh::type_iterator type_it = mesh.firstType(dim, ghost_type); Mesh::type_iterator type_end = mesh.lastType(dim, ghost_type); for (; type_it != type_end; ++type_it) { assembleStiffnessMatrixInterface(type, *type_it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Computes the reinforcement stiffness matrix (Gomes & Awruch, 2001) * \f[ * \mathbf{K}_e = \sum_{i=1}^R{A_i\int_{S_i}{\mathbf{B}^T * \mathbf{C}_i^T \mathbf{D}_{s, i} \mathbf{C}_i \mathbf{B}\,\mathrm{d}s}} * \f] */ template void MaterialReinforcement::assembleStiffnessMatrixInterface( const ElementType & interface_type, const ElementType & background_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt voigt_size = VoigtHelper::size; FEEngine & interface_engine = emodel.getFEEngine("EmbeddedInterfaceFEEngine"); Array & elem_filter = this->element_filter(interface_type, ghost_type); Array & grad_u = gradu_embedded(interface_type, ghost_type); UInt nb_element = elem_filter.size(); UInt nodes_per_background_e = Mesh::getNbNodesPerElement(background_type); UInt nb_quadrature_points = interface_engine.getNbIntegrationPoints(interface_type, ghost_type); UInt back_dof = dim * nodes_per_background_e; UInt integrant_size = back_dof; grad_u.resize(nb_quadrature_points * nb_element); Array tangent_moduli(nb_element * nb_quadrature_points, 1, "interface_tangent_moduli"); this->computeTangentModuli(interface_type, tangent_moduli, ghost_type); Array & shapesd = (*shape_derivatives(interface_type, ghost_type))( background_type, ghost_type); Array integrant(nb_element * nb_quadrature_points, integrant_size * integrant_size, "B^t*C^t*D*C*B"); /// Temporary matrices for integrant product Matrix Bvoigt(voigt_size, back_dof); Matrix DCB(1, back_dof); Matrix CtDCB(voigt_size, back_dof); Array::scalar_iterator D_it = tangent_moduli.begin(); Array::scalar_iterator D_end = tangent_moduli.end(); Array::matrix_iterator C_it = directing_cosines(interface_type, ghost_type).begin(1, voigt_size); Array::matrix_iterator B_it = shapesd.begin(dim, nodes_per_background_e); Array::matrix_iterator integrant_it = integrant.begin(integrant_size, integrant_size); for (; D_it != D_end; ++D_it, ++C_it, ++B_it, ++integrant_it) { Real & D = *D_it; Matrix & C = *C_it; Matrix & B = *B_it; Matrix & BtCtDCB = *integrant_it; VoigtHelper::transferBMatrixToSymVoigtBMatrix(B, Bvoigt, nodes_per_background_e); DCB.mul(C, Bvoigt); DCB *= D * area; CtDCB.mul(C, DCB); BtCtDCB.mul(Bvoigt, CtDCB); } tangent_moduli.resize(0); Array K_interface(nb_element, integrant_size * integrant_size, "K_interface"); interface_engine.integrate(integrant, K_interface, integrant_size * integrant_size, interface_type, ghost_type, elem_filter); integrant.resize(0); // Mauro: Here K_interface contains the local stiffness matrices, // directing_cosines contains the information about the orientation // of the reinforcements, any rotation of the local stiffness matrix // can be done here auto & filter = getBackgroundFilter(interface_type, background_type, ghost_type); emodel.getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", K_interface, background_type, ghost_type, _symmetric, filter); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template Real MaterialReinforcement::getEnergy(const std::string & id) { AKANTU_DEBUG_IN(); if (id == "potential") { Real epot = 0.; this->computePotentialEnergyByElements(); Mesh::type_iterator it = this->element_filter.firstType( this->spatial_dimension), end = this->element_filter.lastType( this->spatial_dimension); for (; it != end; ++it) { FEEngine & interface_engine = emodel.getFEEngine("EmbeddedInterfaceFEEngine"); epot += interface_engine.integrate( this->potential_energy(*it, _not_ghost), *it, _not_ghost, this->element_filter(*it, _not_ghost)); epot *= area; } return epot; } AKANTU_DEBUG_OUT(); return 0; } /* -------------------------------------------------------------------------- */ template void MaterialReinforcement::computeGradU( const ElementType & interface_type, GhostType ghost_type) { // Looping over background types for (auto && bg_type : background_filter(interface_type, ghost_type)->elementTypes(dim)) { const UInt nodes_per_background_e = Mesh::getNbNodesPerElement(bg_type); const UInt voigt_size = VoigtHelper::size; auto & bg_shapesd = (*shape_derivatives(interface_type, ghost_type))(bg_type, ghost_type); auto & filter = getBackgroundFilter(interface_type, bg_type, ghost_type); Array disp_per_element(0, dim * nodes_per_background_e, "disp_elem"); FEEngine::extractNodalToElementField( emodel.getMesh(), emodel.getDisplacement(), disp_per_element, bg_type, ghost_type, filter); Matrix concrete_du(dim, dim); Matrix epsilon(dim, dim); Vector evoigt(voigt_size); for (auto && tuple : zip(make_view(disp_per_element, dim, nodes_per_background_e), make_view(bg_shapesd, dim, nodes_per_background_e), this->gradu(interface_type, ghost_type), make_view(this->directing_cosines(interface_type, ghost_type), voigt_size))) { auto & u = std::get<0>(tuple); auto & B = std::get<1>(tuple); auto & du = std::get<2>(tuple); auto & C = std::get<3>(tuple); concrete_du.mul(u, B); auto epsilon = 0.5 * (concrete_du + concrete_du.transpose()); strainTensorToVoigtVector(epsilon, evoigt); du = C.dot(evoigt); } } } /* -------------------------------------------------------------------------- */ /** * The structure of the directing cosines matrix is : * \f{eqnarray*}{ * C_{1,\cdot} & = & (l^2, m^2, n^2, mn, ln, lm) \\ * C_{i,j} & = & 0 * \f} * * with : * \f[ * (l, m, n) = \frac{1}{\|\frac{\mathrm{d}\vec{r}(s)}{\mathrm{d}s}\|} \cdot * \frac{\mathrm{d}\vec{r}(s)}{\mathrm{d}s} * \f] */ template inline void MaterialReinforcement::computeDirectingCosinesOnQuad( const Matrix & nodes, Matrix & cosines) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(nodes.cols() == 2, "Higher order reinforcement elements not implemented"); const Vector a = nodes(0), b = nodes(1); Vector delta = b - a; Real sq_length = 0.; for (UInt i = 0; i < dim; i++) { sq_length += delta(i) * delta(i); } if (dim == 2) { cosines(0, 0) = delta(0) * delta(0); // l^2 cosines(0, 1) = delta(1) * delta(1); // m^2 cosines(0, 2) = delta(0) * delta(1); // lm } else if (dim == 3) { cosines(0, 0) = delta(0) * delta(0); // l^2 cosines(0, 1) = delta(1) * delta(1); // m^2 cosines(0, 2) = delta(2) * delta(2); // n^2 cosines(0, 3) = delta(1) * delta(2); // mn cosines(0, 4) = delta(0) * delta(2); // ln cosines(0, 5) = delta(0) * delta(1); // lm } cosines /= sq_length; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void MaterialReinforcement::stressTensorToVoigtVector( const Matrix & tensor, Vector & vector) { AKANTU_DEBUG_IN(); for (UInt i = 0; i < dim; i++) { vector(i) = tensor(i, i); } if (dim == 2) { vector(2) = tensor(0, 1); } else if (dim == 3) { vector(3) = tensor(1, 2); vector(4) = tensor(0, 2); vector(5) = tensor(0, 1); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void MaterialReinforcement::strainTensorToVoigtVector( const Matrix & tensor, Vector & vector) { AKANTU_DEBUG_IN(); for (UInt i = 0; i < dim; i++) { vector(i) = tensor(i, i); } if (dim == 2) { vector(2) = 2 * tensor(0, 1); } else if (dim == 3) { vector(3) = 2 * tensor(1, 2); vector(4) = 2 * tensor(0, 2); vector(5) = 2 * tensor(0, 1); } AKANTU_DEBUG_OUT(); } } // akantu diff --git a/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model_prestress.cc b/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model_prestress.cc index 6ed24517d..116108dc1 100644 --- a/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model_prestress.cc +++ b/test/test_model/test_solid_mechanics_model/test_embedded_interface/test_embedded_interface_model_prestress.cc @@ -1,225 +1,224 @@ /** * @file test_embedded_interface_model_prestress.cc * * @author Lucas Frerot * * @date creation: Tue Apr 28 2015 * @date last modification: Thu Oct 15 2015 * * @brief Embedded model test for prestressing (bases on stress norm) * * @section LICENSE * * Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory * (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ #include "aka_common.hh" #include "embedded_interface_model.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; #define YG 0.483644859 #define I_eq 0.012488874 #define A_eq (1e-2 + 1. / 7. * 1.) /* -------------------------------------------------------------------------- */ struct StressSolution : public BC::Neumann::FromHigherDim { Real M; Real I; Real yg; Real pre_stress; StressSolution(UInt dim, Real M, Real I, Real yg = 0, Real pre_stress = 0) : BC::Neumann::FromHigherDim(Matrix(dim, dim)), M(M), I(I), yg(yg), pre_stress(pre_stress) {} virtual ~StressSolution() {} void operator()(const IntegrationPoint & /*quad_point*/, Vector & dual, const Vector & coord, const Vector & normals) const { UInt dim = coord.size(); if (dim < 2) AKANTU_DEBUG_ERROR("Solution not valid for 1D"); Matrix stress(dim, dim); stress.clear(); stress(0, 0) = this->stress(coord(1)); dual.mul(stress, normals); } inline Real stress(Real height) const { return -M / I * (height - yg) + pre_stress; } inline Real neutral_axis() const { return -I * pre_stress / M + yg; } }; /* -------------------------------------------------------------------------- */ int main (int argc, char * argv[]) { initialize("prestress.dat", argc, argv); debug::setDebugLevel(dblError); Math::setTolerance(1e-6); const UInt dim = 2; /* -------------------------------------------------------------------------- */ Mesh mesh(dim); mesh.read("embedded_mesh_prestress.msh"); // mesh.createGroupsFromMeshData("physical_names"); Mesh reinforcement_mesh(dim, "reinforcement_mesh"); try { reinforcement_mesh.read("embedded_mesh_prestress_reinforcement.msh"); } catch(debug::Exception & e) {} // reinforcement_mesh.createGroupsFromMeshData("physical_names"); EmbeddedInterfaceModel model(mesh, reinforcement_mesh, dim); model.initFull(EmbeddedInterfaceModelOptions(_static)); /* -------------------------------------------------------------------------- */ /* Computation of analytical residual */ /* -------------------------------------------------------------------------- */ /* * q = 1000 N/m * L = 20 m * a = 1 m */ - Real steel_area = - model.getMaterial("reinforcement").get("area"); - Real pre_stress = 1e6; + Real steel_area = model.getMaterial("reinforcement").get("area"); + Real pre_stress = model.getMaterial("reinforcement").get("pre_stress"); Real stress_norm = 0.; StressSolution * concrete_stress = nullptr, * steel_stress = nullptr; Real pre_force = pre_stress * steel_area; Real pre_moment = -pre_force * (YG - 0.25); Real neutral_axis = YG - I_eq / A_eq * pre_force / pre_moment; concrete_stress = new StressSolution(dim, pre_moment, 7. * I_eq, YG, -pre_force / (7. * A_eq)); steel_stress = new StressSolution(dim, pre_moment, I_eq, YG, pre_stress - pre_force / A_eq); stress_norm = std::abs(concrete_stress->stress(1)) * (1 - neutral_axis) * 0.5 + std::abs(concrete_stress->stress(0)) * neutral_axis * 0.5 + std::abs(steel_stress->stress(0.25)) * steel_area; model.applyBC(*concrete_stress, "XBlocked"); auto end_node = *mesh.getElementGroup("EndNode").getNodeGroup().begin(); Vector end_node_force = model.getForce().begin(dim)[end_node]; end_node_force(0) += steel_stress->stress(0.25) * steel_area; Array analytical_residual(mesh.getNbNodes(), dim, "analytical_residual"); analytical_residual.copy(model.getForce()); model.getForce().clear(); delete concrete_stress; delete steel_stress; /* -------------------------------------------------------------------------- */ model.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "XBlocked"); model.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "YBlocked"); try { model.solveStep(); } catch (debug::Exception e) { std::cerr << e.what() << std::endl; return EXIT_FAILURE; } /* -------------------------------------------------------------------------- */ /* Computation of FEM residual norm */ /* -------------------------------------------------------------------------- */ ElementGroup & xblocked = mesh.getElementGroup("XBlocked"); NodeGroup & boundary_nodes = xblocked.getNodeGroup(); NodeGroup::const_node_iterator nodes_it = boundary_nodes.begin(), nodes_end = boundary_nodes.end(); model.assembleInternalForces(); Array residual(mesh.getNbNodes(), dim, "my_residual"); residual.copy(model.getInternalForce()); residual -= model.getForce(); Array::vector_iterator com_res = residual.begin(dim); Array::vector_iterator position = mesh.getNodes().begin(dim); Real res_sum = 0.; UInt lower_node = -1; UInt upper_node = -1; Real lower_dist = 1; Real upper_dist = 1; for (; nodes_it != nodes_end ; ++nodes_it) { UInt node_number = *nodes_it; const Vector res = com_res[node_number]; const Vector pos = position[node_number]; if (!Math::are_float_equal(pos(1), 0.25)) { if ((std::abs(pos(1) - 0.25) < lower_dist) && (pos(1) < 0.25)) { lower_dist = std::abs(pos(1) - 0.25); lower_node = node_number; } if ((std::abs(pos(1) - 0.25) < upper_dist) && (pos(1) > 0.25)) { upper_dist = std::abs(pos(1) - 0.25); upper_node = node_number; } } for (UInt i = 0 ; i < dim ; i++) { if (!Math::are_float_equal(pos(1), 0.25)) { res_sum += std::abs(res(i)); } } } const Vector upper_res = com_res[upper_node], lower_res = com_res[lower_node]; const Vector end_node_res = com_res[end_node]; Vector delta = upper_res - lower_res; delta *= lower_dist / (upper_dist + lower_dist); Vector concrete_residual = lower_res + delta; Vector steel_residual = end_node_res - concrete_residual; for (UInt i = 0 ; i < dim ; i++) { res_sum += std::abs(concrete_residual(i)); res_sum += std::abs(steel_residual(i)); } Real relative_error = std::abs(res_sum - stress_norm) / stress_norm; if (relative_error > 1e-3) { std::cerr << "Relative error = " << relative_error << std::endl; return EXIT_FAILURE; } finalize(); return 0; }