diff --git a/src/fe_engine/element_class_structural.hh b/src/fe_engine/element_class_structural.hh index c0f40e9d5..f12f582ad 100644 --- a/src/fe_engine/element_class_structural.hh +++ b/src/fe_engine/element_class_structural.hh @@ -1,201 +1,204 @@ /** * @file element_class_structural.hh * * @author Fabian Barras * @author Nicolas Richart * @author Damien Spielmann * * @date creation: Thu Feb 21 2013 * @date last modification: Thu Oct 22 2015 * * @brief Specialization of the element classes for structural elements * * @section LICENSE * * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "element_class.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_ELEMENT_CLASS_STRUCTURAL_HH__ #define __AKANTU_ELEMENT_CLASS_STRUCTURAL_HH__ namespace akantu { /// Macro to generate the InterpolationProperty structures for different /// interpolation types #define AKANTU_DEFINE_STRUCTURAL_INTERPOLATION_TYPE_PROPERTY( \ itp_type, itp_geom_type, ndof, nb_stress) \ template <> struct InterpolationProperty { \ static const InterpolationKind kind{_itk_structural}; \ static const UInt nb_nodes_per_element{ \ InterpolationProperty::nb_nodes_per_element}; \ static const InterpolationType itp_geometry_type{itp_geom_type}; \ static const UInt natural_space_dimension{ \ InterpolationProperty::natural_space_dimension}; \ static const UInt nb_degree_of_freedom{ndof}; \ static const UInt nb_stress_components{nb_stress}; \ } /* -------------------------------------------------------------------------- */ template class InterpolationElement { public: using interpolation_property = InterpolationProperty; /// compute the shape values for a given set of points in natural coordinates static inline void computeShapes(const Matrix & natural_coord, Tensor3 & N) { for (UInt i = 0; i < natural_coord.cols(); ++i) { Matrix n_t = N(i); computeShapes(natural_coord(i), n_t); } } /// compute the shape values for a given point in natural coordinates static inline void computeShapes(const Vector & natural_coord, Matrix & N); /// compute shape derivatives (input is dxds) for a set of points static inline void computeShapeDerivatives(const Tensor3 & /*Js*/, const Tensor3 & /*DNDSs*/, Tensor3 & /*Bs*/) { AKANTU_DEBUG_TO_IMPLEMENT(); } /** * compute @f$ B_{ij} = \frac{\partial N_j}{\partial S_i} @f$ the variation of * shape functions along with variation of natural coordinates on a given set * of points in natural coordinates */ static inline void computeDNDS(const Matrix & natural_coord, Tensor3 & dnds) { for (UInt i = 0; i < natural_coord.cols(); ++i) { Matrix dnds_t = dnds(i); computeDNDS(natural_coord(i), dnds_t); } } /** * compute @f$ B_{ij} = \frac{\partial N_j}{\partial S_i} @f$ the variation of * shape functions along with * variation of natural coordinates on a given point in natural * coordinates */ static inline void computeDNDS(const Vector & natural_coord, Matrix & dnds); public: static AKANTU_GET_MACRO_NOT_CONST( NbNodesPerInterpolationElement, interpolation_property::nb_nodes_per_element, UInt); static AKANTU_GET_MACRO_NOT_CONST( ShapeSize, (interpolation_property::nb_nodes_per_element * interpolation_property::nb_degree_of_freedom * interpolation_property::nb_degree_of_freedom), UInt); static AKANTU_GET_MACRO_NOT_CONST( ShapeDerivativesSize, (interpolation_property::nb_nodes_per_element * interpolation_property::nb_degree_of_freedom * interpolation_property::nb_stress_components), UInt); static AKANTU_GET_MACRO_NOT_CONST( NaturalSpaceDimension, interpolation_property::natural_space_dimension, UInt); static AKANTU_GET_MACRO_NOT_CONST( NbDegreeOfFreedom, interpolation_property::nb_degree_of_freedom, UInt); static AKANTU_GET_MACRO_NOT_CONST( NbStressComponents, interpolation_property::nb_stress_components, UInt); }; /// Macro to generate the element class structures for different structural /// element types /* -------------------------------------------------------------------------- */ #define AKANTU_DEFINE_STRUCTURAL_ELEMENT_CLASS_PROPERTY( \ elem_type, geom_type, interp_type, parent_el_type, elem_kind, sp, \ gauss_int_type, min_int_order) \ template <> struct ElementClassProperty { \ static const GeometricalType geometrical_type{geom_type}; \ static const InterpolationType interpolation_type{interp_type}; \ static const ElementType parent_element_type{parent_el_type}; \ static const ElementKind element_kind{elem_kind}; \ static const UInt spatial_dimension{sp}; \ static const GaussIntegrationType gauss_integration_type{gauss_int_type}; \ static const UInt polynomial_degree{min_int_order}; \ } /* -------------------------------------------------------------------------- */ /* ElementClass for structural elements */ /* -------------------------------------------------------------------------- */ template class ElementClass : public GeometricalElement< ElementClassProperty::geometrical_type>, public InterpolationElement< ElementClassProperty::interpolation_type> { protected: using geometrical_element = GeometricalElement::geometrical_type>; using interpolation_element = InterpolationElement< ElementClassProperty::interpolation_type>; using parent_element = ElementClass::parent_element_type>; public: /// compute jacobian (or integration variable change factor) for a given point static inline void computeJMat(const Tensor3 & /*DNDSs*/, const Matrix & /*Xs*/, Tensor3 & /*Js*/) { AKANTU_DEBUG_TO_IMPLEMENT(); } static inline void computeJacobian(const Matrix & /*natural_coords*/, const Matrix & /*node_coords*/, Vector & /*jacobians*/) { AKANTU_DEBUG_TO_IMPLEMENT(); } + static inline void computeRotation(const Matrix & node_coords, + Matrix & rotation); + public: static AKANTU_GET_MACRO_NOT_CONST(Kind, _ek_structural, ElementKind); static AKANTU_GET_MACRO_NOT_CONST(P1ElementType, _not_defined, ElementType); static AKANTU_GET_MACRO_NOT_CONST(FacetType, _not_defined, ElementType); static ElementType getFacetType(__attribute__((unused)) UInt t = 0) { return _not_defined; } static AKANTU_GET_MACRO_NOT_CONST( SpatialDimension, ElementClassProperty::spatial_dimension, UInt); static ElementType * getFacetTypeInternal() { return NULL; } }; } // namespace akantu /* -------------------------------------------------------------------------- */ #include "element_classes/element_class_hermite_inline_impl.cc" /* keep order */ #include "element_classes/element_class_bernoulli_beam_inline_impl.cc" #include "element_classes/element_class_kirchhoff_shell_inline_impl.cc" /* -------------------------------------------------------------------------- */ #endif /* __AKANTU_ELEMENT_CLASS_STRUCTURAL_HH__ */ diff --git a/src/fe_engine/element_classes/element_class_bernoulli_beam_inline_impl.cc b/src/fe_engine/element_classes/element_class_bernoulli_beam_inline_impl.cc index 6ae762fb3..1887abfaa 100644 --- a/src/fe_engine/element_classes/element_class_bernoulli_beam_inline_impl.cc +++ b/src/fe_engine/element_classes/element_class_bernoulli_beam_inline_impl.cc @@ -1,152 +1,165 @@ /** * @file element_class_bernoulli_beam_inline_impl.cc * * @author Fabian Barras * * @date creation: Fri Jul 15 2011 * @date last modification: Sun Oct 19 2014 * * @brief Specialization of the element_class class for the type _bernoulli_beam_2 * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * * @section DESCRIPTION * * @verbatim --x-----q1----|----q2-----x---> x -1 0 1 @endverbatim * */ /* -------------------------------------------------------------------------- */ #include "aka_static_if.hh" #include "element_class_structural.hh" //#include "aka_element_classes_info.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_ELEMENT_CLASS_BERNOULLI_BEAM_INLINE_IMPL_CC__ #define __AKANTU_ELEMENT_CLASS_BERNOULLI_BEAM_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ AKANTU_DEFINE_STRUCTURAL_INTERPOLATION_TYPE_PROPERTY(_itp_bernoulli_beam_2, _itp_lagrange_segment_2, 3, 2); AKANTU_DEFINE_STRUCTURAL_INTERPOLATION_TYPE_PROPERTY(_itp_bernoulli_beam_3, _itp_lagrange_segment_2, 6, 4); AKANTU_DEFINE_STRUCTURAL_ELEMENT_CLASS_PROPERTY(_bernoulli_beam_2, _gt_segment_2, _itp_bernoulli_beam_2, _segment_2, _ek_structural, 2, _git_segment, 5); AKANTU_DEFINE_STRUCTURAL_ELEMENT_CLASS_PROPERTY(_bernoulli_beam_3, _gt_segment_2, _itp_bernoulli_beam_3, _segment_2, _ek_structural, 3, _git_segment, 5); /* -------------------------------------------------------------------------- */ template <> inline void InterpolationElement<_itp_bernoulli_beam_2, _itk_structural>::computeShapes( const Vector & natural_coords, Matrix & N) { Vector L(2); InterpolationElement<_itp_lagrange_segment_2, _itk_lagrangian>::computeShapes( natural_coords, L); Matrix H(2, 4); InterpolationElement<_itp_hermite_2, _itk_structural>::computeShapes( natural_coords, H); // clang-format off // u1 v1 t1 u2 v2 t2 N = {{L(0), 0 , 0 , L(1), 0 , 0 }, // u {0 , H(0, 0), H(0, 1), 0 , H(0, 2), H(0, 3)}, // v {0 , H(1, 0), H(1, 1), 0 , H(1, 2), H(1, 3)}}; // theta // clang-format on } template <> inline void InterpolationElement<_itp_bernoulli_beam_3, _itk_structural>::computeShapes( const Vector & natural_coords, Matrix & N) { Vector L(2); InterpolationElement<_itp_lagrange_segment_2, _itk_lagrangian>::computeShapes( natural_coords, L); Matrix H(2, 4); InterpolationElement<_itp_hermite_2, _itk_structural>::computeShapes( natural_coords, H); // clang-format off // u1 v1 w1 x1 y1 z1 u2 v2 w2 x2 y2 z2 N = {{L(0), 0 , 0 , 0 , 0 , 0 , L(1), 0 , 0 , 0 , 0 , 0 }, // u {0 , H(0, 0), 0 , 0 , H(0, 1), 0 , 0 , H(0, 2), 0 , 0 , H(0, 3), 0 }, // v {0 , 0 , H(0, 0), 0 , 0 , H(0, 1), 0 , 0 , H(0, 2), 0 , 0 , H(0, 3)}, // w {0 , 0 , 0 , L(0), 0 , 0 , 0 , 0 , 0 , L(1), 0 , 0 }, // thetax {0 , H(1, 0), 0 , 0 , H(1, 1), 0 , 0 , H(1, 2), 0 , 0 , H(1, 3), 0 }, // thetay {0 , 0 , H(1, 0), 0 , 0 , H(1, 1), 0 , 0 , H(1, 2), 0 , 0 , H(1, 3)}}; // thetaz // clang-format on } /* -------------------------------------------------------------------------- */ template <> inline void InterpolationElement<_itp_bernoulli_beam_2, _itk_structural>::computeDNDS( const Vector & natural_coords, Matrix & B) { Matrix L(1, 2); InterpolationElement<_itp_lagrange_segment_2, _itk_lagrangian>::computeDNDS( natural_coords, L); Matrix H(1, 4); InterpolationElement<_itp_hermite_2, _itk_structural>::computeDNDS( natural_coords, H); // clang-format off // u1 v1 t1 u2 v2 t2 B = {{L(0, 0), 0 , 0 , L(0, 1), 0 , 0 }, // epsilon {0 , H(0, 0), H(0, 1), 0 , H(0, 2), H(0, 3)}}; // chi (curv.) // clang-format on } template <> inline void InterpolationElement<_itp_bernoulli_beam_3, _itk_structural>::computeDNDS( const Vector & natural_coords, Matrix & B) { Matrix L(1, 2); InterpolationElement<_itp_lagrange_segment_2, _itk_lagrangian>::computeDNDS( natural_coords, L); Matrix H(1, 4); InterpolationElement<_itp_hermite_2, _itk_structural>::computeDNDS( natural_coords, H); // clang-format off // u1 v1 w1 x1 y1 z1 u2 v2 w2 x2 y2 z2 B = {{L(0, 0), 0 , 0 , 0 , 0 , 0 , L(0, 1), 0 , 0 , 0 , 0 , 0 }, // eps {0 , H(0, 0), 0 , 0 , H(0, 1), 0 , 0 , H(0, 2), 0 , 0 , H(0, 3), 0 }, // chix {0 , 0 , H(0, 0), 0 , 0 , H(0, 1), 0 , 0 , H(0, 2), 0 , 0 , H(0, 3)}, // chiy {0 , 0 , 0 , L(0, 0), 0 , 0 , 0 , 0 , 0 , L(0, 1), 0 , 0 }}; // chiz // clang-format on } +template <> +inline void ElementClass<_bernoulli_beam_2>::computeRotation( + const Matrix & node_coords, Matrix & rotation) { + auto X1 = node_coords(0); + auto X2 = node_coords(1); + auto vec = Vector(X2) - Vector(X1); + auto L = vec.norm(); + auto c = vec(0) / L; + auto s = vec(1) / L; + + rotation = {{c, -s}, {s, c}}; +} + } // namespace akantu #endif /* __AKANTU_ELEMENT_CLASS_BERNOULLI_BEAM_INLINE_IMPL_CC__ */ diff --git a/src/fe_engine/shape_structural.hh b/src/fe_engine/shape_structural.hh index d3d14d331..dca244784 100644 --- a/src/fe_engine/shape_structural.hh +++ b/src/fe_engine/shape_structural.hh @@ -1,114 +1,119 @@ /** * @file shape_structural.hh * * @author Fabian Barras * @author Nicolas Richart * * @date creation: Tue Feb 15 2011 * @date last modification: Thu Oct 22 2015 * * @brief shape class for element with different set of shapes functions * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "shape_functions.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SHAPE_STRUCTURAL_HH__ #define __AKANTU_SHAPE_STRUCTURAL_HH__ namespace akantu { template class ShapeStructural : public ShapeFunctions { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef ElementTypeMap **> ElementTypeMapMultiReal; ShapeStructural(Mesh & mesh, const ID & id = "shape_structural", const MemoryID & memory_id = 0); virtual ~ShapeStructural(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// compute shape functions on given integration points template void computeShapesOnIntegrationPoints( const Array &, const Matrix & integration_points, Array & shapes, const GhostType & ghost_type, const Array & filter_elements = empty_filter) const; /// initialization function for structural elements inline void initShapeFunctions(const Array & nodes, const Matrix & integration_points, const ElementType & type, const GhostType & ghost_type); /// pre compute all shapes on the element integration points from natural /// coordinates template void precomputeShapesOnIntegrationPoints(const Array & nodes, const GhostType & ghost_type); /// pre compute all shapes on the element integration points from natural /// coordinates template void precomputeShapeDerivativesOnIntegrationPoints(const Array & nodes, const GhostType & ghost_type); + /// pre compute all rotation matrices + template + void precomputeRotationMatrices(const Array & nodes, + const GhostType & ghost_type); + /// interpolate nodal values on the integration points template void interpolateOnIntegrationPoints( const Array & u, Array & uq, UInt nb_degree_of_freedom, const GhostType & ghost_type = _not_ghost, const Array & filter_elements = empty_filter) const; /// compute the gradient of u on the integration points template void gradientOnIntegrationPoints( const Array & u, Array & nablauq, UInt nb_degree_of_freedom, const GhostType & ghost_type = _not_ghost, const Array & filter_elements = empty_filter) const; /// multiply a field by shape functions template void fieldTimesShapes(__attribute__((unused)) const Array & field, __attribute__((unused)) Array & fiedl_times_shapes, __attribute__((unused)) const GhostType & ghost_type) const { AKANTU_DEBUG_TO_IMPLEMENT(); } protected: ElementTypeMapArray rotation_matrices; }; } // namespace akantu #include "shape_structural_inline_impl.cc" #endif /* __AKANTU_SHAPE_STRUCTURAL_HH__ */ diff --git a/src/fe_engine/shape_structural_inline_impl.cc b/src/fe_engine/shape_structural_inline_impl.cc index a83064837..273003c7e 100644 --- a/src/fe_engine/shape_structural_inline_impl.cc +++ b/src/fe_engine/shape_structural_inline_impl.cc @@ -1,308 +1,331 @@ /** * @file shape_structural_inline_impl.cc * * @author Fabian Barras * @author Nicolas Richart * * @date creation: Mon Dec 13 2010 * @date last modification: Thu Oct 15 2015 * * @brief ShapeStructural inline implementation * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "mesh_iterators.hh" #include "shape_structural.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SHAPE_STRUCTURAL_INLINE_IMPL_CC__ #define __AKANTU_SHAPE_STRUCTURAL_INLINE_IMPL_CC__ namespace akantu { template inline void ShapeStructural::initShapeFunctions( const Array & /* unused */, const Matrix & /* unused */, const ElementType & /* unused */, const GhostType & /* unused */) { AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ #define INIT_SHAPE_FUNCTIONS(type) \ setIntegrationPointsByType(integration_points, ghost_type); \ precomputeShapesOnIntegrationPoints(nodes, ghost_type); \ precomputeShapeDerivativesOnIntegrationPoints(nodes, ghost_type); template <> inline void ShapeStructural<_ek_structural>::initShapeFunctions( const Array & nodes, const Matrix & integration_points, const ElementType & type, const GhostType & ghost_type) { AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(INIT_SHAPE_FUNCTIONS); } #undef INIT_SHAPE_FUNCTIONS /* -------------------------------------------------------------------------- */ template <> template void ShapeStructural<_ek_structural>::computeShapesOnIntegrationPoints( const Array & /*nodes*/, const Matrix & integration_points, Array & shapes, const GhostType & ghost_type, const Array & filter_elements) const { /// \TODO this code differs from ShapeLagrangeBase only in the size of N UInt nb_points = integration_points.cols(); UInt nb_element = mesh.getConnectivity(type, ghost_type).size(); shapes.resize(nb_element * nb_points); UInt ndof = ElementClass::getNbDegreeOfFreedom(); #if !defined(AKANTU_NDEBUG) UInt size_of_shapes = ElementClass::getShapeSize(); AKANTU_DEBUG_ASSERT(shapes.getNbComponent() == size_of_shapes, "The shapes array does not have the correct " << "number of component"); #endif auto shapes_it = shapes.begin_reinterpret( ElementClass::getNbNodesPerInterpolationElement(), ndof, nb_points, nb_element); auto shapes_begin = shapes_it; if (filter_elements != empty_filter) { nb_element = filter_elements.size(); } for (UInt elem = 0; elem < nb_element; ++elem) { if (filter_elements != empty_filter) shapes_it = shapes_begin + filter_elements(elem); Tensor3 & N = *shapes_it; ElementClass::computeShapes(integration_points, N); if (filter_elements == empty_filter) ++shapes_it; } } /* -------------------------------------------------------------------------- */ template <> template void ShapeStructural<_ek_structural>::precomputeShapesOnIntegrationPoints( const Array & /*nodes*/, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); const auto & natural_coords = integration_points(type, ghost_type); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_points = integration_points(type, ghost_type).cols(); auto nb_element = mesh.getNbElement(type, ghost_type); auto nb_dof = ElementClass::getNbDegreeOfFreedom(); auto itp_type = FEEngine::getInterpolationType(type); if (not shapes.exists(itp_type, ghost_type)) { auto size_of_shapes = this->getShapeSize(type); this->shapes.alloc(0, size_of_shapes, itp_type, ghost_type); } auto & shapes_ = this->shapes(itp_type, ghost_type); shapes_.resize(nb_element * nb_points); auto shapes_it = shapes_.begin_reinterpret( nb_dof, nb_dof * nb_nodes_per_element, nb_points, nb_element); for (UInt elem = 0; elem < nb_element; ++elem, ++shapes_it) { auto & N = *shapes_it; ElementClass::computeShapes(natural_coords, N); } AKANTU_DEBUG_OUT(); } // namespace akantu /* -------------------------------------------------------------------------- */ template template void ShapeStructural::precomputeShapeDerivativesOnIntegrationPoints( const Array & nodes, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); const auto & natural_coords = integration_points(type, ghost_type); auto spatial_dimension = mesh.getSpatialDimension(); auto natural_spatial_dimension = ElementClass::getNaturalSpaceDimension(); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_points = natural_coords.cols(); auto nb_dof = ElementClass::getNbDegreeOfFreedom(); auto nb_element = mesh.getNbElement(type, ghost_type); auto nb_stress_components = ElementClass::getNbStressComponents(); auto itp_type = FEEngine::getInterpolationType(type); if (not this->shapes_derivatives.exists(itp_type, ghost_type)) { auto size_of_shapesd = this->getShapeDerivativesSize(type); this->shapes_derivatives.alloc(0, size_of_shapesd, itp_type, ghost_type); } auto rot_matrices = this->rotation_matrices(type, ghost_type); Array x_el(0, spatial_dimension * nb_nodes_per_element); FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type); auto & shapesd = this->shapes_derivatives(itp_type, ghost_type); shapesd.resize(nb_element * nb_points); for (auto && tuple : zip(make_view(x_el, spatial_dimension, nb_nodes_per_element), make_view(shapesd, nb_stress_components, nb_nodes_per_element * nb_dof, nb_points), make_view(rot_matrices, spatial_dimension, spatial_dimension))) { // compute shape derivatives auto & X = std::get<0>(tuple); auto & B = std::get<1>(tuple); auto & R = std::get<2>(tuple); // Rotate to local basis auto x = R * X; Matrix node_coords(natural_spatial_dimension, nb_nodes_per_element); // Extract relevant first lines for (UInt j = 0; j < node_coords.rows(); ++j) for (UInt i = 0; i < node_coords.cols(); ++i) node_coords(i, j) = x(i, j); Tensor3 dnds(B.size(0), B.size(1), B.size(2)); ElementClass::computeDNDS(natural_coords, dnds); Tensor3 J(node_coords.rows(), natural_coords.rows(), natural_coords.cols()); ElementClass::computeJMat(dnds, node_coords, J); ElementClass::computeShapeDerivatives(J, dnds, B); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void ShapeStructural::interpolateOnIntegrationPoints( const Array & in_u, Array & out_uq, UInt nb_dof, const GhostType & ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(out_uq.getNbComponent() == nb_dof, "The output array shape is not correct"); auto itp_type = FEEngine::getInterpolationType(type); const auto & shapes_ = shapes(itp_type, ghost_type); auto nb_element = mesh.getNbElement(type, ghost_type); auto nb_nodes_per_element = ElementClass::getNbNodesPerElement(); auto nb_quad_points_per_element = integration_points(type, ghost_type).cols(); Array u_el(0, nb_nodes_per_element * nb_dof); FEEngine::extractNodalToElementField(mesh, in_u, u_el, ghost_type, filter_elements); auto nb_quad_points = nb_quad_points_per_element * u_el.size(); out_uq.resize(nb_quad_points); auto out_it = out_uq.begin_reinterpret(nb_dof, 1, nb_quad_points_per_element, u_el.size()); auto shapes_it = shapes_.begin_reinterpret(nb_dof, nb_dof * nb_nodes_per_element, nb_quad_points_per_element, nb_element); auto u_it = u_el.begin_reinterpret(nb_dof * nb_nodes_per_element, 1, nb_quad_points_per_element, u_el.size()); for_each_elements(nb_element, filter_elements, [&](auto && el) { auto & uq = *out_it; const auto & u = *u_it; auto N = Tensor3(shapes_it[el]); for (auto && q : arange(uq.size(2))) { auto uq_q = Matrix(uq(q)); auto u_q = Matrix(u(q)); auto N_q = Matrix(N(q)); uq_q.mul(N, u); } ++out_it; ++u_it; }); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void ShapeStructural::gradientOnIntegrationPoints( const Array & in_u, Array & out_nablauq, UInt nb_dof, const GhostType & ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); auto itp_type = FEEngine::getInterpolationType(type); const auto & shapesd = shapes_derivatives(itp_type, ghost_type); auto nb_element = mesh.getNbElement(type, ghost_type); auto element_dimension = ElementClass::getSpatialDimension(); auto nb_quad_points_per_element = integration_points(type, ghost_type).cols(); auto nb_nodes_per_element = ElementClass::getNbNodesPerElement(); Array u_el(0, nb_nodes_per_element * nb_dof); FEEngine::extractNodalToElementField(mesh, in_u, u_el, ghost_type, filter_elements); auto nb_quad_points = nb_quad_points_per_element * u_el.size(); out_nablauq.resize(nb_quad_points); auto out_it = out_nablauq.begin_reinterpret( element_dimension, 1, nb_quad_points_per_element, u_el.size()); auto shapesd_it = shapesd.begin_reinterpret( element_dimension, nb_dof * nb_nodes_per_element, nb_quad_points_per_element, nb_element); auto u_it = u_el.begin_reinterpret(nb_dof * nb_nodes_per_element, 1, nb_quad_points_per_element, u_el.size()); for_each_elements(nb_element, filter_elements, [&](auto && el) { auto & nablau = *out_it; const auto & u = *u_it; auto B = Tensor3(shapesd_it[el]); for (auto && q : arange(nablau.size(2))) { auto nablau_q = Matrix(uq(q)); auto u_q = Matrix(u(q)); auto B_q = Matrix(N(q)); nablau_q.mul(B, u); } ++out_it; ++u_it; }); AKANTU_DEBUG_OUT(); } +template +template +void ShapeStructural::precomputeRotationMatrices( + const Array & nodes, const GhostType & ghost_type) { + auto element_dimension = ElementClass::getSpatialDimension(); + auto nb_nodes_per_element = ElementClass::getNbNodesPerElement(); + + Array x_el(0, nb_nodes_per_element * element_dimension); + FEEngine::extractNodalToElementField(mesh, nodes, x_el, ghost_type); + + rotation_matrices(type).resize(element_dimension * element_dimension * + x_el.size()); + + for (auto && tuple : + zip(make_view(rotation_matrices(type), element_dimension, + element_dimension), + make_view(x_el, element_dimension, nb_nodes_per_element))) { + auto & R = std::get<0>(tuple); + auto & X = std::get<1>(tuple); + ElementClass::computeRotation(X, R); + } +} + } // namespace akantu #endif /* __AKANTU_SHAPE_STRUCTURAL_INLINE_IMPL_CC__ */ diff --git a/src/model/structural_mechanics/structural_elements/structural_element_bernoulli_beam_2.hh b/src/model/structural_mechanics/structural_elements/structural_element_bernoulli_beam_2.hh index 0dd653cf6..5086ec2fa 100644 --- a/src/model/structural_mechanics/structural_elements/structural_element_bernoulli_beam_2.hh +++ b/src/model/structural_mechanics/structural_elements/structural_element_bernoulli_beam_2.hh @@ -1,143 +1,134 @@ /** * @file structural_element_bernoulli_beam_2.hh * * @author Fabian Barras * @author Sébastien Hartmann * @author Nicolas Richart * * @date creation Tue Sep 19 2017 * * @brief Specific functions for bernoulli beam 2d * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "structural_mechanics_model.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_STRUCTURAL_ELEMENT_BERNOULLI_BEAM_2_HH__ #define __AKANTU_STRUCTURAL_ELEMENT_BERNOULLI_BEAM_2_HH__ namespace akantu { /* -------------------------------------------------------------------------- */ template <> inline void StructuralMechanicsModel::assembleMass<_bernoulli_beam_2>() { AKANTU_DEBUG_IN(); -#if 1 constexpr ElementType type = _bernoulli_beam_2; auto & fem = getFEEngineClass(); auto nb_element = mesh.getNbElement(type); auto nb_nodes_per_element = mesh.getNbNodesPerElement(type); auto nb_quadrature_points = fem.getNbIntegrationPoints(type); auto nb_fields_to_interpolate = ElementClass::getNbStressComponents(); auto nt_n_field_size = nb_degree_of_freedom * nb_nodes_per_element; Array n(nb_element * nb_quadrature_points, nb_fields_to_interpolate * nt_n_field_size, "N"); - //const auto & N_star = fem.getShapes(type); - - // for - Array * rho_field = new Array(nb_element * nb_quadrature_points, 1, "Rho"); rho_field->clear(); computeRho(*rho_field, type, _not_ghost); +#if 0 bool sign = true; - //for (auto && ghost_type : ghost_types) { - // fem.computeShapesMatrix(type, nb_degree_of_freedom, nb_nodes_per_element, - // n, - // 0, 0, 0, sign, ghost_type); // Ni ui -> u - // fem.computeShapesMatrix(type, nb_degree_of_freedom, nb_nodes_per_element, - // n, - // 1, 1, 1, sign, ghost_type); // Mi vi -> v - // fem.computeShapesMatrix(type, nb_degree_of_freedom, nb_nodes_per_element, - // n, - // 2, 2, 1, sign, ghost_type); // Li Theta_i -> v - // fem.assembleFieldMatrix(*rho_field, nb_degree_of_freedom, *mass_matrix, - // n, - // rotation_matrix, type, ghost_type); - //} - - // delete n; + for (auto && ghost_type : ghost_types) { + fem.computeShapesMatrix(type, nb_degree_of_freedom, nb_nodes_per_element, n, + 0, 0, 0, sign, ghost_type); // Ni ui -> u + fem.computeShapesMatrix(type, nb_degree_of_freedom, nb_nodes_per_element, n, + 1, 1, 1, sign, ghost_type); // Mi vi -> v + fem.computeShapesMatrix(type, nb_degree_of_freedom, nb_nodes_per_element, n, + 2, 2, 1, sign, ghost_type); // Li Theta_i -> v + fem.assembleFieldMatrix(*rho_field, nb_degree_of_freedom, *mass_matrix, n, + rotation_matrix, type, ghost_type); + } +#endif + delete rho_field; AKANTU_DEBUG_OUT(); -#endif } /* -------------------------------------------------------------------------- */ template <> void StructuralMechanicsModel::computeRotationMatrix<_bernoulli_beam_2>( Array & rotations) { auto type = _bernoulli_beam_2; auto nb_element = mesh.getNbElement(type); auto connec_it = mesh.getConnectivity(type).begin(2); auto nodes_it = mesh.getNodes().begin(spatial_dimension); auto R_it = rotations.begin(nb_degree_of_freedom, nb_degree_of_freedom); for (UInt e = 0; e < nb_element; ++e, ++R_it, ++connec_it) { auto & R = *R_it; auto & connec = *connec_it; Vector x2 = nodes_it[connec(1)]; // X2 Vector x1 = nodes_it[connec(0)]; // X1 auto le = x1.distance(x2); auto c = (x2(0) - x1(0)) / le; auto s = (x2(1) - x1(1)) / le; /// Definition of the rotation matrix R = {{c, s, 0.}, {-s, c, 0.}, {0., 0., 1.}}; } } /* -------------------------------------------------------------------------- */ template <> void StructuralMechanicsModel::computeTangentModuli<_bernoulli_beam_2>( Array & tangent_moduli) { //auto nb_element = getFEEngine().getMesh().getNbElement(_bernoulli_beam_2); auto nb_quadrature_points = getFEEngine().getNbIntegrationPoints(_bernoulli_beam_2); auto tangent_size = 2; auto D_it = tangent_moduli.begin(tangent_size, tangent_size); auto el_mat = element_material(_bernoulli_beam_2, _not_ghost).begin(); for (auto & mat : element_material(_bernoulli_beam_2, _not_ghost)) { auto E = materials[mat].E; auto A = materials[mat].A; auto I = materials[mat].I; for (UInt q = 0; q < nb_quadrature_points; ++q, ++D_it) { auto & D = *D_it; D(0, 0) = E * A; D(1, 1) = E * I; } } } } // namespace akantu #endif /* __AKANTU_STRUCTURAL_ELEMENT_BERNOULLI_BEAM_2_HH__ */ diff --git a/src/model/structural_mechanics/structural_mechanics_model_inline_impl.cc b/src/model/structural_mechanics/structural_mechanics_model_inline_impl.cc index d8a5f99c1..ec27c7ede 100644 --- a/src/model/structural_mechanics/structural_mechanics_model_inline_impl.cc +++ b/src/model/structural_mechanics/structural_mechanics_model_inline_impl.cc @@ -1,363 +1,361 @@ /** * @file structural_mechanics_model_inline_impl.cc * * @author Fabian Barras * @author Sébastien Hartmann * @author Nicolas Richart * @author Damien Spielmann * * @date creation: Fri Jul 15 2011 * @date last modification: Thu Oct 15 2015 * * @brief Implementation of inline functions of StructuralMechanicsModel * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "structural_mechanics_model.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_STRUCTURAL_MECHANICS_MODEL_INLINE_IMPL_CC__ #define __AKANTU_STRUCTURAL_MECHANICS_MODEL_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ template void StructuralMechanicsModel::computeTangentModuli( Array & /*tangent_moduli*/) { AKANTU_DEBUG_TO_IMPLEMENT(); } } // namespace akantu #include "structural_element_bernoulli_beam_2.hh" #include "structural_element_bernoulli_beam_3.hh" #include "structural_element_kirchhoff_shell.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template void StructuralMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); auto nb_element = getFEEngine().getMesh().getNbElement(type); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type); auto tangent_size = ElementClass::getNbStressComponents(); auto tangent_moduli = std::make_unique>( nb_element * nb_quadrature_points, tangent_size * tangent_size, "tangent_stiffness_matrix"); computeTangentModuli(*tangent_moduli); /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_d_b_size = nb_degree_of_freedom * nb_nodes_per_element; auto bt_d_b = std::make_unique>( nb_element * nb_quadrature_points, bt_d_b_size * bt_d_b_size, "B^t*D*B"); const auto & b = getFEEngine().getShapesDerivatives(type); - Matrix Bt_D(bt_d_b_size, tangent_size); - - auto B = b.begin(tangent_size, bt_d_b_size); - auto D = tangent_moduli->begin(tangent_size, tangent_size); - auto Bt_D_B = bt_d_b->begin(bt_d_b_size, bt_d_b_size); - auto T = rotation_matrix(type).begin(bt_d_b_size, bt_d_b_size); - - for (UInt e = 0; e < nb_element; ++e, ++T) { - for (UInt q = 0; q < nb_quadrature_points; ++q, ++B, ++D, ++Bt_D_B) { - auto BT = *B * *T; - Bt_D.mul(BT, *D); - Bt_D_B->mul(Bt_D, BT); - } + Matrix BtD(bt_d_b_size, tangent_size); + + for (auto && tuple : + zip(make_view(b, tangent_size, bt_d_b_size), + make_view(*tangent_moduli, tangent_size, tangent_size), + make_view(*bt_d_b, bt_d_b_size, bt_d_b_size))) { + auto & B = std::get<0>(tuple); + auto & D = std::get<1>(tuple); + auto & BtDB = std::get<2>(tuple); + BtD.mul(B, D); + BtDB.template mul(BtD, B); } /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ auto int_bt_d_b = std::make_unique>( nb_element, bt_d_b_size * bt_d_b_size, "int_B^t*D*B"); getFEEngine().integrate(*bt_d_b, *int_bt_d_b, bt_d_b_size * bt_d_b_size, type); getDOFManager().assembleElementalMatricesToMatrix( "K", "displacement", *int_bt_d_b, type, _not_ghost, _symmetric); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void StructuralMechanicsModel::computeStressOnQuad() { AKANTU_DEBUG_IN(); Array & sigma = stress(type, _not_ghost); auto nb_element = mesh.getNbElement(type); auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type); auto tangent_size = ElementClass::getNbStressComponents(); auto tangent_moduli = std::make_unique>( nb_element * nb_quadrature_points, tangent_size * tangent_size, "tangent_stiffness_matrix"); computeTangentModuli(*tangent_moduli); /// compute DB auto d_b_size = nb_degree_of_freedom * nb_nodes_per_element; auto d_b = std::make_unique>(nb_element * nb_quadrature_points, d_b_size * tangent_size, "D*B"); const auto & b = getFEEngine().getShapesDerivatives(type); auto B = b.begin(tangent_size, d_b_size); auto D = tangent_moduli->begin(tangent_size, tangent_size); auto D_B = d_b->begin(tangent_size, d_b_size); for (UInt e = 0; e < nb_element; ++e) { for (UInt q = 0; q < nb_quadrature_points; ++q, ++B, ++D, ++D_B) { D_B->template mul(*D, *B); } } /// compute DBu D_B = d_b->begin(tangent_size, d_b_size); auto DBu = sigma.begin(tangent_size); Vector ul(d_b_size); Array u_el(0, d_b_size); FEEngine::extractNodalToElementField(mesh, *displacement_rotation, u_el, type); auto ug = u_el.begin(d_b_size); auto T = rotation_matrix(type).begin(d_b_size, d_b_size); for (UInt e = 0; e < nb_element; ++e, ++T, ++ug) { ul.mul(*T, *ug); for (UInt q = 0; q < nb_quadrature_points; ++q, ++D_B, ++DBu) { DBu->template mul(*D_B, ul); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void StructuralMechanicsModel::computeForcesByLocalTractionArray( const Array & tractions) { AKANTU_DEBUG_IN(); #if 0 UInt nb_element = getFEEngine().getMesh().getNbElement(type); UInt nb_nodes_per_element = getFEEngine().getMesh().getNbNodesPerElement(type); UInt nb_quad = getFEEngine().getNbIntegrationPoints(type); // check dimension match AKANTU_DEBUG_ASSERT( Mesh::getSpatialDimension(type) == getFEEngine().getElementDimension(), "element type dimension does not match the dimension of boundaries : " << getFEEngine().getElementDimension() << " != " << Mesh::getSpatialDimension(type)); // check size of the vector AKANTU_DEBUG_ASSERT( tractions.size() == nb_quad * nb_element, "the size of the vector should be the total number of quadrature points"); // check number of components AKANTU_DEBUG_ASSERT(tractions.getNbComponent() == nb_degree_of_freedom, "the number of components should be the spatial " "dimension of the problem"); Array Nvoigt(nb_element * nb_quad, nb_degree_of_freedom * nb_degree_of_freedom * nb_nodes_per_element); transferNMatrixToSymVoigtNMatrix(Nvoigt); Array::const_matrix_iterator N_it = Nvoigt.begin( nb_degree_of_freedom, nb_degree_of_freedom * nb_nodes_per_element); Array::const_matrix_iterator T_it = rotation_matrix(type).begin(nb_degree_of_freedom * nb_nodes_per_element, nb_degree_of_freedom * nb_nodes_per_element); Array::const_vector_iterator te_it = tractions.begin(nb_degree_of_freedom); Array funct(nb_element * nb_quad, nb_degree_of_freedom * nb_nodes_per_element, 0.); Array::iterator> Fe_it = funct.begin(nb_degree_of_freedom * nb_nodes_per_element); Vector fe(nb_degree_of_freedom * nb_nodes_per_element); for (UInt e = 0; e < nb_element; ++e, ++T_it) { const Matrix & T = *T_it; for (UInt q = 0; q < nb_quad; ++q, ++N_it, ++te_it, ++Fe_it) { const Matrix & N = *N_it; const Vector & te = *te_it; Vector & Fe = *Fe_it; // compute N^t tl fe.mul(N, te); // turn N^t tl back in the global referential Fe.mul(T, fe); } } // allocate the vector that will contain the integrated values std::stringstream name; name << id << type << ":integral_boundary"; Array int_funct(nb_element, nb_degree_of_freedom * nb_nodes_per_element, name.str()); // do the integration getFEEngine().integrate(funct, int_funct, nb_degree_of_freedom * nb_nodes_per_element, type); // assemble the result into force vector getFEEngine().assembleArray(int_funct, *force_momentum, dof_synchronizer->getLocalDOFEquationNumbers(), nb_degree_of_freedom, type); #endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void StructuralMechanicsModel::computeForcesByGlobalTractionArray( const Array & traction_global) { AKANTU_DEBUG_IN(); #if 0 UInt nb_element = mesh.getNbElement(type); UInt nb_quad = getFEEngine().getNbIntegrationPoints(type); UInt nb_nodes_per_element = getFEEngine().getMesh().getNbNodesPerElement(type); std::stringstream name; name << id << ":structuralmechanics:imposed_linear_load"; Array traction_local(nb_element * nb_quad, nb_degree_of_freedom, name.str()); Array::const_matrix_iterator T_it = rotation_matrix(type).begin(nb_degree_of_freedom * nb_nodes_per_element, nb_degree_of_freedom * nb_nodes_per_element); Array::const_iterator> Te_it = traction_global.begin(nb_degree_of_freedom); Array::iterator> te_it = traction_local.begin(nb_degree_of_freedom); Matrix R(nb_degree_of_freedom, nb_degree_of_freedom); for (UInt e = 0; e < nb_element; ++e, ++T_it) { const Matrix & T = *T_it; for (UInt i = 0; i < nb_degree_of_freedom; ++i) for (UInt j = 0; j < nb_degree_of_freedom; ++j) R(i, j) = T(i, j); for (UInt q = 0; q < nb_quad; ++q, ++Te_it, ++te_it) { const Vector & Te = *Te_it; Vector & te = *te_it; // turn the traction in the local referential te.mul(R, Te); } } computeForcesByLocalTractionArray(traction_local); #endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * @param myf pointer to a function that fills a vector/tensor with respect to * passed coordinates */ #if 0 template inline void StructuralMechanicsModel::computeForcesFromFunction( BoundaryFunction myf, BoundaryFunctionType function_type) { /** function type is ** _bft_forces : linear load is given ** _bft_stress : stress function is given -> Not already done for this kind *of model */ std::stringstream name; name << id << ":structuralmechanics:imposed_linear_load"; Array lin_load(0, nb_degree_of_freedom, name.str()); name.clear(); UInt offset = nb_degree_of_freedom; // prepare the loop over element types UInt nb_quad = getFEEngine().getNbIntegrationPoints(type); UInt nb_element = getFEEngine().getMesh().getNbElement(type); name.clear(); name << id << ":structuralmechanics:quad_coords"; Array quad_coords(nb_element * nb_quad, spatial_dimension, "quad_coords"); getFEEngineClass() .getShapeFunctions() .interpolateOnIntegrationPoints(getFEEngine().getMesh().getNodes(), quad_coords, spatial_dimension); getFEEngineClass() .getShapeFunctions() .interpolateOnIntegrationPoints( getFEEngine().getMesh().getNodes(), quad_coords, spatial_dimension, _not_ghost, empty_filter, true, 0, 1, 1); if (spatial_dimension == 3) getFEEngineClass() .getShapeFunctions() .interpolateOnIntegrationPoints( getFEEngine().getMesh().getNodes(), quad_coords, spatial_dimension, _not_ghost, empty_filter, true, 0, 2, 2); lin_load.resize(nb_element * nb_quad); Real * imposed_val = lin_load.storage(); /// sigma/load on each quadrature points Real * qcoord = quad_coords.storage(); for (UInt el = 0; el < nb_element; ++el) { for (UInt q = 0; q < nb_quad; ++q) { myf(qcoord, imposed_val, NULL, 0); imposed_val += offset; qcoord += spatial_dimension; } } switch (function_type) { case _bft_traction_local: computeForcesByLocalTractionArray(lin_load); break; case _bft_traction: computeForcesByGlobalTractionArray(lin_load); break; default: break; } } #endif } // namespace akantu #endif /* __AKANTU_STRUCTURAL_MECHANICS_MODEL_INLINE_IMPL_CC__ */