diff --git a/src/fe_engine/element_classes/element_class_kirchhoff_shell_inline_impl.cc b/src/fe_engine/element_classes/element_class_kirchhoff_shell_inline_impl.cc index 6818c8ddc..4769bf237 100644 --- a/src/fe_engine/element_classes/element_class_kirchhoff_shell_inline_impl.cc +++ b/src/fe_engine/element_classes/element_class_kirchhoff_shell_inline_impl.cc @@ -1,211 +1,211 @@ /** * @file element_class_kirchhoff_shell_inline_impl.cc * * @author Damien Spielmann * * @date creation: Fri Jul 04 2014 * @date last modification: Sun Oct 19 2014 * * @brief Element class Kirchhoff Shell * * @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_structural.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_ELEMENT_CLASS_KIRCHHOFF_SHELL_INLINE_IMPL_CC__ #define __AKANTU_ELEMENT_CLASS_KIRCHHOFF_SHELL_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ AKANTU_DEFINE_STRUCTURAL_INTERPOLATION_TYPE_PROPERTY( _itp_discrete_kirchhoff_triangle_18, _itp_lagrange_triangle_3, 6, 6, 21); AKANTU_DEFINE_STRUCTURAL_ELEMENT_CLASS_PROPERTY( _discrete_kirchhoff_triangle_18, _gt_triangle_3, _itp_discrete_kirchhoff_triangle_18, _triangle_3, _ek_structural, 3, _git_triangle, 2); /* -------------------------------------------------------------------------- */ namespace detail { inline void computeBasisChangeMatrix(Matrix & P, const Matrix & X) { Vector X1 = X(0); Vector X2 = X(1); Vector X3 = X(2); Vector a1 = X2 - X1; Vector a2 = X3 - X1; a1.normalize(); Vector e3 = a1.crossProduct(a2); e3.normalize(); Vector e2 = e3.crossProduct(a1); P(0) = a1; P(1) = e2; P(2) = e3; } } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_discrete_kirchhoff_triangle_18>::computeRotationMatrix( Matrix & R, const Matrix & X, const Vector &) { auto dim = X.rows(); Matrix P(dim, dim); detail::computeBasisChangeMatrix(P, X); R.clear(); for (UInt i = 0; i < dim; ++i) for (UInt j = 0; j < dim; ++j) R(i + dim, j + dim) = R(i, j) = P(j, i); } /* -------------------------------------------------------------------------- */ template <> inline void InterpolationElement<_itp_discrete_kirchhoff_triangle_18>::computeShapes( const Vector & /*natural_coords*/, const Matrix & /*real_coord*/, Matrix & /*N*/) {} /* -------------------------------------------------------------------------- */ template <> inline void InterpolationElement<_itp_discrete_kirchhoff_triangle_18>::computeDNDS( const Vector & natural_coords, const Matrix & real_coordinates, Matrix & B) { auto dim = real_coordinates.cols(); Matrix P(dim, dim); detail::computeBasisChangeMatrix(P, real_coordinates); auto X = P * real_coordinates; Vector X1 = X(0); Vector X2 = X(1); Vector X3 = X(2); std::array, 3> A = {X2 - X1, X3 - X2, X1 - X3}; std::array L, C, S; // Setting all last coordinates to 0 std::for_each(A.begin(), A.end(), [](auto & a) { a(2) = 0; }); // Computing lengths std::transform(A.begin(), A.end(), L.begin(), [](Vector & a) { return a.norm(); }); // Computing cosines std::transform(A.begin(), A.end(), L.begin(), C.begin(), [](Vector & a, Real & l) { return a(0) / l; }); // Computing sines std::transform(A.begin(), A.end(), L.begin(), S.begin(), [](Vector & a, Real & l) { return a(1) / l; }); // Natural coordinates Real xi = natural_coords(0); Real eta = natural_coords(1); // Derivative of quadratic interpolation functions Matrix dP = {{4 * (1 - 2 * xi - eta), 4 * eta, -4 * eta}, {-4 * xi, 4 * xi, 4 * (1 - xi - 2 * eta)}}; Matrix dNx1 = { {3. / 2 * (dP(0, 0) * C[0] / L[0] - dP(0, 2) * C[2] / L[2]), 3. / 2 * (dP(0, 1) * C[1] / L[1] - dP(0, 0) * C[0] / L[0]), 3. / 2 * (dP(0, 2) * C[2] / L[2] - dP(0, 1) * C[1] / L[1])}, {3. / 2 * (dP(1, 0) * C[0] / L[0] - dP(1, 2) * C[2] / L[2]), 3. / 2 * (dP(1, 1) * C[1] / L[1] - dP(1, 0) * C[0] / L[0]), 3. / 2 * (dP(1, 2) * C[2] / L[2] - dP(1, 1) * C[1] / L[1])}}; Matrix dNx2 = { // clang-format off {-1 - 3. / 4 * (dP(0, 0) * C[0] * C[0] + dP(0, 2) * C[2] * C[2]), 1 - 3. / 4 * (dP(0, 1) * C[1] * C[1] + dP(0, 0) * C[0] * C[0]), -3. / 4 * (dP(0, 2) * C[2] * C[2] + dP(0, 1) * C[1] * C[1])}, {-1 - 3. / 4 * (dP(1, 0) * C[0] * C[0] + dP(1, 2) * C[2] * C[2]), -3. / 4 * (dP(1, 1) * C[1] * C[1] + dP(1, 0) * C[0] * C[0]), 1 - 3. / 4 * (dP(1, 2) * C[2] * C[2] + dP(1, 1) * C[1] * C[1])}}; // clang-format on Matrix dNx3 = { {-3. / 4 * (dP(0, 0) * C[0] * S[0] + dP(0, 2) * C[2] * S[2]), -3. / 4 * (dP(0, 1) * C[1] * S[1] + dP(0, 0) * C[0] * S[0]), -3. / 4 * (dP(0, 2) * C[2] * S[2] + dP(0, 1) * C[1] * S[1])}, {-3. / 4 * (dP(1, 0) * C[0] * S[0] + dP(1, 2) * C[2] * S[2]), -3. / 4 * (dP(1, 1) * C[1] * S[1] + dP(1, 0) * C[0] * S[0]), -3. / 4 * (dP(1, 2) * C[2] * S[2] + dP(1, 1) * C[1] * S[1])}}; Matrix dNy1 = { {3. / 2 * (dP(0, 0) * S[0] / L[0] - dP(0, 2) * S[2] / L[2]), 3. / 2 * (dP(0, 1) * S[1] / L[1] - dP(0, 0) * S[0] / L[0]), 3. / 2 * (dP(0, 2) * S[2] / L[2] - dP(0, 1) * S[1] / L[1])}, {3. / 2 * (dP(1, 0) * S[0] / L[0] - dP(1, 2) * S[2] / L[2]), 3. / 2 * (dP(1, 1) * S[1] / L[1] - dP(1, 0) * S[0] / L[0]), 3. / 2 * (dP(1, 2) * S[2] / L[2] - dP(1, 1) * S[1] / L[1])}}; Matrix dNy2 = dNx3; Matrix dNy3 = { // clang-format off {-1 - 3. / 4 * (dP(0, 0) * S[0] * S[0] + dP(0, 2) * S[2] * S[2]), 1 - 3. / 4 * (dP(0, 1) * S[1] * S[1] + dP(0, 0) * S[0] * S[0]), -3. / 4 * (dP(0, 2) * S[2] * S[2] + dP(0, 1) * S[1] * S[1])}, {-1 - 3. / 4 * (dP(1, 0) * S[0] * S[0] + dP(1, 2) * S[2] * S[2]), -3. / 4 * (dP(1, 1) * S[1] * S[1] + dP(1, 0) * S[0] * S[0]), 1 - 3. / 4 * (dP(1, 2) * S[2] * S[2] + dP(1, 1) * S[1] * S[1])}}; // clang-format on // Derivative of linear (membrane mode) functions Matrix dNm(2, 3); InterpolationElement<_itp_lagrange_triangle_3, _itk_lagrangian>::computeDNDS( natural_coords, dNm); UInt i = 0; for (const Matrix & mat : {dNm, dNx1, dNx2, dNx3, dNy1, dNy2, dNy3}) { B.block(mat, 0, i); i += mat.cols(); } } /* -------------------------------------------------------------------------- */ template <> inline void InterpolationElement<_itp_discrete_kirchhoff_triangle_18, _itk_structural>::arrangeInVoigt(const Matrix & dnds, Matrix & B) { Matrix dNm(2, 3), dNx1(2, 3), dNx2(2, 3), dNx3(2, 3), dNy1(2, 3), dNy2(2, 3), dNy3(2, 3); UInt i = 0; for (Matrix * mat : {&dNm, &dNx1, &dNx2, &dNx3, &dNy1, &dNy2, &dNy3}) { *mat = dnds.block(0, i, 2, 3); i += mat->cols(); } // clang-format off B = { // Membrane shape functions; TODO use the triangle_3 shapes - {dNm(0, 0), 0, 0, 0, 0, 0, dNm(0, 1), 0, 0, 0, 0, 0, dNm(0, 2), 0, 0, 0, 0, 0, }, - {0, dNm(1, 0), 0, 0, 0, 0, 0, dNm(1, 1), 0, 0, 0, 0, 0, dNm(1, 2), 0, 0, 0, 0, }, - {dNm(1, 0), dNm(0, 0), 0, 0, 0, 0, dNm(1, 1), dNm(0, 1), 0, 0, 0, 0, dNm(1, 2), dNm(0, 2), 0, 0, 0, 0, }, + {dNm(0, 0), 0, 0, 0, 0, 0, dNm(0, 1), 0, 0, 0, 0, 0, dNm(0, 2), 0, 0, 0, 0, 0}, + {0, dNm(1, 0), 0, 0, 0, 0, 0, dNm(1, 1), 0, 0, 0, 0, 0, dNm(1, 2), 0, 0, 0, 0}, + {dNm(1, 0), dNm(0, 0), 0, 0, 0, 0, dNm(1, 1), dNm(0, 1), 0, 0, 0, 0, dNm(1, 2), dNm(0, 2), 0, 0, 0, 0}, // Bending shape functions {0, 0, dNx1(0, 0), -dNx3(0, 0), dNx2(0, 0), 0, 0, 0, dNx1(0, 1), -dNx3(0, 1), dNx2(0, 1), 0, 0, 0, dNx1(0, 2), -dNx3(0, 2), dNx2(0, 2), 0}, {0, 0, dNy1(1, 0), -dNy3(1, 0), dNy2(1, 0), 0, 0, 0, dNy1(1, 1), -dNy3(1, 1), dNy2(1, 1), 0, 0, 0, dNy1(1, 2), -dNy3(1, 2), dNy2(1, 2), 0}, {0, 0, dNx1(1, 0) + dNy1(0, 0), -dNx3(1, 0) - dNy3(0, 0), dNx2(1, 0) + dNy2(1, 0), 0, 0, 0, dNx1(1, 1) + dNy1(0, 1), -dNx3(1, 1) - dNy3(0, 0), dNx2(1, 1) + dNy2(1, 1), 0, 0, 0, dNx1(1, 2) + dNy1(0, 2), -dNx3(1, 2) - dNy3(0, 2), dNx2(1, 2) + dNy2(1, 2), 0}}; // clang-format on } } // namespace akantu #endif /* __AKANTU_ELEMENT_CLASS_KIRCHHOFF_SHELL_INLINE_IMPL_CC__ */ diff --git a/src/fe_engine/shape_structural_inline_impl.cc b/src/fe_engine/shape_structural_inline_impl.cc index 5792b0d34..207ab5bea 100644 --- a/src/fe_engine/shape_structural_inline_impl.cc +++ b/src/fe_engine/shape_structural_inline_impl.cc @@ -1,389 +1,431 @@ /** * @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 { namespace { /// Extract nodal coordinates per elements template std::unique_ptr> getNodesPerElement(const Mesh & mesh, const Array & nodes, const GhostType & ghost_type) { const auto dim = ElementClass::getSpatialDimension(); const auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); auto nodes_per_element = std::make_unique>(0, dim * nb_nodes_per_element); FEEngine::extractNodalToElementField(mesh, nodes, *nodes_per_element, type, ghost_type); return nodes_per_element; } } 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); \ precomputeRotationMatrices(nodes, 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 { 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(); } auto nodes_per_element = getNodesPerElement(mesh, nodes, ghost_type); auto nodes_it = nodes_per_element->begin(mesh.getSpatialDimension(), Mesh::getNbNodesPerElement(type)); auto nodes_begin = nodes_it; for (UInt elem = 0; elem < nb_element; ++elem) { if (filter_elements != empty_filter) { shapes_it = shapes_begin + filter_elements(elem); nodes_it = nodes_begin + filter_elements(elem); } Tensor3 & N = *shapes_it; auto & real_coord = *nodes_it; ElementClass::computeShapes(integration_points, real_coord, N); if (filter_elements == empty_filter) ++shapes_it; } } /* -------------------------------------------------------------------------- */ template template void ShapeStructural::precomputeRotationMatrices( const Array & nodes, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); const auto spatial_dimension = mesh.getSpatialDimension(); const auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const auto nb_element = mesh.getNbElement(type, ghost_type); const auto nb_dof = ElementClass::getNbDegreeOfFreedom(); if (not this->rotation_matrices.exists(type, ghost_type)) { this->rotation_matrices.alloc(0, nb_dof * nb_dof, type, ghost_type); } auto & rot_matrices = this->rotation_matrices(type, ghost_type); rot_matrices.resize(nb_element); Array x_el(0, spatial_dimension * nb_nodes_per_element); FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type); bool has_extra_normal = mesh.hasData("extra_normal", type, ghost_type); Array::const_vector_iterator extra_normal; if (has_extra_normal) extra_normal = mesh.getData("extra_normal", type, ghost_type) .begin(spatial_dimension); for (auto && tuple : zip(make_view(x_el, spatial_dimension, nb_nodes_per_element), make_view(rot_matrices, nb_dof, nb_dof))) { // compute shape derivatives auto & X = std::get<0>(tuple); auto & R = std::get<1>(tuple); if (has_extra_normal) { ElementClass::computeRotationMatrix(R, X, *extra_normal); ++extra_normal; } else { ElementClass::computeRotationMatrix( R, X, Vector(spatial_dimension)); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void ShapeStructural::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(); const auto dim = ElementClass::getSpatialDimension(); 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 nodes_per_element = getNodesPerElement(mesh, nodes, ghost_type); for (auto && tuple : zip(make_view(shapes_, nb_dof, nb_dof * nb_nodes_per_element, nb_points), make_view(*nodes_per_element, dim, nb_nodes_per_element))) { auto & N = std::get<0>(tuple); auto & real_coord = std::get<1>(tuple); ElementClass::computeShapes(natural_coords, real_coord, N); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void ShapeStructural::precomputeShapeDerivativesOnIntegrationPoints( const Array & nodes, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); const auto & natural_coords = integration_points(type, ghost_type); const auto spatial_dimension = mesh.getSpatialDimension(); const auto natural_spatial_dimension = ElementClass::getNaturalSpaceDimension(); const auto nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const auto nb_points = natural_coords.cols(); const auto nb_dof = ElementClass::getNbDegreeOfFreedom(); const auto nb_element = mesh.getNbElement(type, ghost_type); const 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, nb_dof, nb_dof))) { // compute shape derivatives auto & X = std::get<0>(tuple); auto & B = std::get<1>(tuple); auto & RDOFs = std::get<2>(tuple); // /!\ This is a temporary variable with no specified size for columns ! // It is up to the element to resize Tensor3 dnds(natural_spatial_dimension, ElementClass::interpolation_property::dnds_columns, B.size(2)); ElementClass::computeDNDS(natural_coords, X, dnds); Tensor3 J(natural_spatial_dimension, natural_coords.rows(), natural_coords.cols()); // Computing the coordinates of the element in the natural space auto R = RDOFs.block(0, 0, spatial_dimension, spatial_dimension); Matrix T(B.size(1), B.size(1)); T.block(RDOFs, 0, 0); T.block(RDOFs, RDOFs.rows(), RDOFs.rows()); // Rotate to local basis auto x = (R * X).block(0, 0, natural_spatial_dimension, nb_nodes_per_element); ElementClass::computeJMat(natural_coords, x, J); ElementClass::computeShapeDerivatives(J, dnds, T, 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, type, 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_element(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_q, u_q); } ++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, type, 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_element(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(nablau(q)); auto u_q = Matrix(u(q)); auto B_q = Matrix(B(q)); nablau_q.mul(B_q, u_q); } ++out_it; ++u_it; }); AKANTU_DEBUG_OUT(); } +/* -------------------------------------------------------------------------- */ +template <> +template +void ShapeStructural<_ek_structural>::computeBtD( + const Array & Ds, Array & BtDs, + GhostType ghost_type, const Array & filter_elements) const { + auto itp_type = ElementClassProperty::interpolation_type; + + auto nb_stress = ElementClass::getNbStressComponents(); + auto nb_dof_per_element = ElementClass::getNbDegreeOfFreedom() * + mesh.getNbNodesPerElement(type); + + const auto & shapes_derivatives = + this->shapes_derivatives(itp_type, ghost_type); + + Array shapes_derivatives_filtered(0, + shapes_derivatives.getNbComponent()); + auto && view = make_view(shapes_derivatives, nb_stress, nb_dof_per_element); + auto B_it = view.begin(); + auto B_end = view.end(); + + if (filter_elements != empty_filter) { + FEEngine::filterElementalData(this->mesh, shapes_derivatives, + shapes_derivatives_filtered, type, ghost_type, + filter_elements); + auto && view = + make_view(shapes_derivatives_filtered, nb_stress, nb_dof_per_element); + B_it = view.begin(); + B_end = view.end(); + } + + for (auto && values : + zip(range(B_it, B_end), + make_view(Ds, nb_stress), + make_view(BtDs, BtDs.getNbComponent()))) { + const auto & B = std::get<0>(values); + const auto & D = std::get<1>(values); + auto & Bt_D = std::get<2>(values); + Bt_D.template mul(B, D); + } +} + } // namespace akantu #endif /* __AKANTU_SHAPE_STRUCTURAL_INLINE_IMPL_CC__ */ diff --git a/src/mesh/element_type_map.cc b/src/mesh/element_type_map.cc index 10ed4c342..ee0fc8cca 100644 --- a/src/mesh/element_type_map.cc +++ b/src/mesh/element_type_map.cc @@ -1,59 +1,71 @@ /** * @file element_type_map.cc * * @author Nicolas Richart * * @date creation Tue Jun 27 2017 * * @brief A Documented file. * * @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 "fe_engine.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ namespace akantu { -FEEngineElementTypeMapArrayInializer::FEEngineElementTypeMapArrayInializer( +FEEngineElementTypeMapArrayInitializer::FEEngineElementTypeMapArrayInitializer( const FEEngine & fe_engine, UInt nb_component, UInt spatial_dimension, const GhostType & ghost_type, const ElementKind & element_kind) - : MeshElementTypeMapArrayInializer( - fe_engine.getMesh(), - nb_component, + : MeshElementTypeMapArrayInitializer( + fe_engine.getMesh(), nb_component, spatial_dimension == UInt(-2) ? fe_engine.getMesh().getSpatialDimension() : spatial_dimension, ghost_type, element_kind, true, false), fe_engine(fe_engine) {} -UInt FEEngineElementTypeMapArrayInializer::size( +FEEngineElementTypeMapArrayInitializer::FEEngineElementTypeMapArrayInitializer( + const FEEngine & fe_engine, + const ElementTypeMapArrayInitializer::CompFunc & nb_component, + UInt spatial_dimension, const GhostType & ghost_type, + const ElementKind & element_kind) + : MeshElementTypeMapArrayInitializer( + fe_engine.getMesh(), nb_component, + spatial_dimension == UInt(-2) + ? fe_engine.getMesh().getSpatialDimension() + : spatial_dimension, + ghost_type, element_kind, true, false), + fe_engine(fe_engine) {} + +UInt FEEngineElementTypeMapArrayInitializer::size( const ElementType & type) const { - return MeshElementTypeMapArrayInializer::size(type) * + return MeshElementTypeMapArrayInitializer::size(type) * fe_engine.getNbIntegrationPoints(type, this->ghost_type); } -FEEngineElementTypeMapArrayInializer::ElementTypesIteratorHelper -FEEngineElementTypeMapArrayInializer::elementTypes() const { +FEEngineElementTypeMapArrayInitializer::ElementTypesIteratorHelper +FEEngineElementTypeMapArrayInitializer::elementTypes() const { return this->fe_engine.elementTypes(spatial_dimension, ghost_type, element_kind); } } // namespace akantu diff --git a/src/mesh/element_type_map.hh b/src/mesh/element_type_map.hh index 2b132a0ee..f3aefd368 100644 --- a/src/mesh/element_type_map.hh +++ b/src/mesh/element_type_map.hh @@ -1,444 +1,449 @@ /** * @file element_type_map.hh * * @author Nicolas Richart * * @date creation: Wed Aug 31 2011 * @date last modification: Fri Oct 02 2015 * * @brief storage class by element type * * @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 "aka_array.hh" #include "aka_memory.hh" #include "aka_named_argument.hh" #include "element.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_ELEMENT_TYPE_MAP_HH__ #define __AKANTU_ELEMENT_TYPE_MAP_HH__ namespace akantu { class FEEngine; } // namespace akantu namespace akantu { namespace { DECLARE_NAMED_ARGUMENT(all_ghost_types); DECLARE_NAMED_ARGUMENT(default_value); DECLARE_NAMED_ARGUMENT(element_kind); DECLARE_NAMED_ARGUMENT(ghost_type); DECLARE_NAMED_ARGUMENT(nb_component); DECLARE_NAMED_ARGUMENT(nb_component_functor); DECLARE_NAMED_ARGUMENT(with_nb_element); DECLARE_NAMED_ARGUMENT(with_nb_nodes_per_element); DECLARE_NAMED_ARGUMENT(spatial_dimension); DECLARE_NAMED_ARGUMENT(do_not_default); } // namespace template class ElementTypeMap; /* -------------------------------------------------------------------------- */ /* ElementTypeMapBase */ /* -------------------------------------------------------------------------- */ /// Common non templated base class for the ElementTypeMap class class ElementTypeMapBase { public: virtual ~ElementTypeMapBase() = default; }; /* -------------------------------------------------------------------------- */ /* ElementTypeMap */ /* -------------------------------------------------------------------------- */ template class ElementTypeMap : public ElementTypeMapBase { public: ElementTypeMap(); ~ElementTypeMap() override; inline static std::string printType(const SupportType & type, const GhostType & ghost_type); /*! Tests whether a type is present in the object * @param type the type to check for * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @return true if the type is present. */ inline bool exists(const SupportType & type, const GhostType & ghost_type = _not_ghost) const; /*! get the stored data corresponding to a type * @param type the type to check for * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @return stored data corresponding to type. */ inline const Stored & operator()(const SupportType & type, const GhostType & ghost_type = _not_ghost) const; /*! get the stored data corresponding to a type * @param type the type to check for * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @return stored data corresponding to type. */ inline Stored & operator()(const SupportType & type, const GhostType & ghost_type = _not_ghost); /*! insert data of a new type (not yet present) into the map. THIS METHOD IS * NOT ARRAY SAFE, when using ElementTypeMapArray, use setArray instead * @param data to insert * @param type type of data (if this type is already present in the map, * an exception is thrown). * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @return stored data corresponding to type. */ template inline Stored & operator()(U && insertee, const SupportType & type, const GhostType & ghost_type = _not_ghost); public: /// print helper virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* Element type Iterator */ /* ------------------------------------------------------------------------ */ /*! iterator allows to iterate over type-data pairs of the map. The interface * expects the SupportType to be ElementType. */ using DataMap = std::map; class type_iterator : private std::iterator { public: using value_type = const SupportType; using pointer = const SupportType *; using reference = const SupportType &; protected: using DataMapIterator = typename ElementTypeMap::DataMap::const_iterator; public: type_iterator(DataMapIterator & list_begin, DataMapIterator & list_end, UInt dim, ElementKind ek); type_iterator(const type_iterator & it); type_iterator() = default; inline reference operator*(); inline reference operator*() const; inline type_iterator & operator++(); type_iterator operator++(int); inline bool operator==(const type_iterator & other) const; inline bool operator!=(const type_iterator & other) const; type_iterator & operator=(const type_iterator & other); private: DataMapIterator list_begin; DataMapIterator list_end; UInt dim; ElementKind kind; }; /// helper class to use in range for constructions class ElementTypesIteratorHelper { public: using Container = ElementTypeMap; using iterator = typename Container::type_iterator; ElementTypesIteratorHelper(const Container & container, UInt dim, GhostType ghost_type, ElementKind kind) : container(std::cref(container)), dim(dim), ghost_type(ghost_type), kind(kind) {} template ElementTypesIteratorHelper(const Container & container, use_named_args_t, pack &&... _pack) : ElementTypesIteratorHelper( container, OPTIONAL_NAMED_ARG(spatial_dimension, _all_dimensions), OPTIONAL_NAMED_ARG(ghost_type, _not_ghost), OPTIONAL_NAMED_ARG(element_kind, _ek_regular)) {} ElementTypesIteratorHelper(const ElementTypesIteratorHelper &) = default; ElementTypesIteratorHelper & operator=(const ElementTypesIteratorHelper &) = default; ElementTypesIteratorHelper & operator=(ElementTypesIteratorHelper &&) = default; iterator begin() { return container.get().firstType(dim, ghost_type, kind); } iterator end() { return container.get().lastType(dim, ghost_type, kind); } private: std::reference_wrapper container; UInt dim; GhostType ghost_type; ElementKind kind; }; private: ElementTypesIteratorHelper elementTypesImpl(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const; template ElementTypesIteratorHelper elementTypesImpl(const use_named_args_t & /*unused*/, pack &&... _pack) const; public: + /*! + * \param _spatial_dimension optional: filter for elements of given spatial + * dimension + * \param _ghost_type optional: filter for a certain ghost_type + * \param _element_kind optional: filter for elements of given kind + */ template std::enable_if_t::value, ElementTypesIteratorHelper> elementTypes(pack &&... _pack) const { return elementTypesImpl(use_named_args, std::forward(_pack)...); } template std::enable_if_t::value, ElementTypesIteratorHelper> elementTypes(pack &&... _pack) const { return elementTypesImpl(std::forward(_pack)...); } public: /*! Get an iterator to the beginning of a subset datamap. This method expects * the SupportType to be ElementType. * @param dim optional: iterate over data of dimension dim (e.g. when * iterating over (surface) facets of a 3D mesh, dim would be 2). * by default, all dimensions are considered. * @param ghost_type optional: by default, the data map for non-ghost * elements is iterated over. * @param kind optional: the kind of element to search for (see * aka_common.hh), by default all kinds are considered * @return an iterator to the first stored data matching the filters * or an iterator to the end of the map if none match*/ inline type_iterator firstType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const; /*! Get an iterator to the end of a subset datamap. This method expects * the SupportType to be ElementType. * @param dim optional: iterate over data of dimension dim (e.g. when * iterating over (surface) facets of a 3D mesh, dim would be 2). * by default, all dimensions are considered. * @param ghost_type optional: by default, the data map for non-ghost * elements is iterated over. * @param kind optional: the kind of element to search for (see * aka_common.hh), by default all kinds are considered * @return an iterator to the last stored data matching the filters * or an iterator to the end of the map if none match */ inline type_iterator lastType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const; protected: /*! Direct access to the underlying data map. for internal use by daughter * classes only * @param ghost_type whether to return the data map or the ghost_data map * @return the raw map */ inline DataMap & getData(GhostType ghost_type); /*! Direct access to the underlying data map. for internal use by daughter * classes only * @param ghost_type whether to return the data map or the ghost_data map * @return the raw map */ inline const DataMap & getData(GhostType ghost_type) const; /* ------------------------------------------------------------------------ */ protected: DataMap data; DataMap ghost_data; }; /* -------------------------------------------------------------------------- */ /* Some typedefs */ /* -------------------------------------------------------------------------- */ template class ElementTypeMapArray : public ElementTypeMap *, SupportType>, public Memory { public: using type = T; using array_type = Array; protected: using parent = ElementTypeMap *, SupportType>; using DataMap = typename parent::DataMap; public: using type_iterator = typename parent::type_iterator; /// standard assigment (copy) operator void operator=(const ElementTypeMapArray &) = delete; ElementTypeMapArray(const ElementTypeMapArray &) = delete; /// explicit copy void copy(const ElementTypeMapArray & other); /*! Constructor * @param id optional: identifier (string) * @param parent_id optional: parent identifier. for organizational purposes * only * @param memory_id optional: choose a specific memory, defaults to memory 0 */ ElementTypeMapArray(const ID & id = "by_element_type_array", const ID & parent_id = "no_parent", const MemoryID & memory_id = 0) : parent(), Memory(parent_id + ":" + id, memory_id), name(id){}; /*! allocate memory for a new array * @param size number of tuples of the new array * @param nb_component tuple size * @param type the type under which the array is indexed in the map * @param ghost_type whether to add the field to the data map or the * ghost_data map * @return a reference to the allocated array */ inline Array & alloc(UInt size, UInt nb_component, const SupportType & type, const GhostType & ghost_type, const T & default_value = T()); /*! allocate memory for a new array in both the data and the ghost_data map * @param size number of tuples of the new array * @param nb_component tuple size * @param type the type under which the array is indexed in the map*/ inline void alloc(UInt size, UInt nb_component, const SupportType & type, const T & default_value = T()); /* get a reference to the array of certain type * @param type data filed under type is returned * @param ghost_type optional: by default the non-ghost map is searched * @return a reference to the array */ inline const Array & operator()(const SupportType & type, const GhostType & ghost_type = _not_ghost) const; /// access the data of an element, this combine the map and array accessor inline const T & operator()(const Element & element, UInt component = 0) const; /// access the data of an element, this combine the map and array accessor inline T & operator()(const Element & element, UInt component = 0); /* get a reference to the array of certain type * @param type data filed under type is returned * @param ghost_type optional: by default the non-ghost map is searched * @return a const reference to the array */ inline Array & operator()(const SupportType & type, const GhostType & ghost_type = _not_ghost); /*! insert data of a new type (not yet present) into the map. * @param type type of data (if this type is already present in the map, * an exception is thrown). * @param ghost_type optional: by default, the data map for non-ghost * elements is searched * @param vect the vector to include into the map * @return stored data corresponding to type. */ inline void setArray(const SupportType & type, const GhostType & ghost_type, const Array & vect); /*! frees all memory related to the data*/ inline void free(); /*! set all values in the ElementTypeMap to zero*/ inline void clear(); /*! set all values in the ElementTypeMap to value */ template inline void set(const ST & value); /*! deletes and reorders entries in the stored arrays * @param new_numbering a ElementTypeMapArray of new indices. UInt(-1) * indicates * deleted entries. */ inline void onElementsRemoved(const ElementTypeMapArray & new_numbering); /// text output helper void printself(std::ostream & stream, int indent = 0) const override; /*! set the id * @param id the new name */ inline void setID(const ID & id) { this->id = id; } ElementTypeMap getNbComponents(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const { ElementTypeMap nb_components; for (auto & type : this->elementTypes(dim, ghost_type, kind)) { UInt nb_comp = (*this)(type, ghost_type).getNbComponent(); nb_components(type, ghost_type) = nb_comp; } return nb_components; } /* ------------------------------------------------------------------------ */ /* more evolved allocators */ /* ------------------------------------------------------------------------ */ public: /// initialize the arrays in accordance to a functor - template - void initialize(const Func & f, const T & default_value, bool do_not_default, - CompFunc && func); + template + void initialize(const Func & f, const T & default_value, bool do_not_default); /// initialize with sizes and number of components in accordance of a mesh /// content template void initialize(const Mesh & mesh, pack &&... _pack); /// initialize with sizes and number of components in accordance of a fe /// engine content (aka integration points) template void initialize(const FEEngine & fe_engine, pack &&... _pack); /* ------------------------------------------------------------------------ */ /* Accesssors */ /* ------------------------------------------------------------------------ */ public: /// get the name of the internal field AKANTU_GET_MACRO(Name, name, ID); /// name of the elment type map: e.g. connectivity, grad_u ID name; }; /// to store data Array by element type using ElementTypeMapReal = ElementTypeMapArray; /// to store data Array by element type using ElementTypeMapInt = ElementTypeMapArray; /// to store data Array by element type using ElementTypeMapUInt = ElementTypeMapArray; /// Map of data of type UInt stored in a mesh using UIntDataMap = std::map *>; using ElementTypeMapUIntDataMap = ElementTypeMap; } // namespace akantu #endif /* __AKANTU_ELEMENT_TYPE_MAP_HH__ */ diff --git a/src/mesh/element_type_map_tmpl.hh b/src/mesh/element_type_map_tmpl.hh index cf1cdc109..23a34e1fe 100644 --- a/src/mesh/element_type_map_tmpl.hh +++ b/src/mesh/element_type_map_tmpl.hh @@ -1,696 +1,747 @@ /** * @file element_type_map_tmpl.hh * * @author Nicolas Richart * * @date creation: Wed Aug 31 2011 * @date last modification: Fri Oct 02 2015 * * @brief implementation of template functions of the ElementTypeMap and * ElementTypeMapArray classes * * @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 "aka_static_if.hh" #include "element_type_map.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ #include "element_type_conversion.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_ELEMENT_TYPE_MAP_TMPL_HH__ #define __AKANTU_ELEMENT_TYPE_MAP_TMPL_HH__ namespace akantu { /* -------------------------------------------------------------------------- */ /* ElementTypeMap */ /* -------------------------------------------------------------------------- */ template inline std::string ElementTypeMap::printType(const SupportType & type, const GhostType & ghost_type) { std::stringstream sstr; sstr << "(" << ghost_type << ":" << type << ")"; return sstr.str(); } /* -------------------------------------------------------------------------- */ template inline bool ElementTypeMap::exists( const SupportType & type, const GhostType & ghost_type) const { return this->getData(ghost_type).find(type) != this->getData(ghost_type).end(); } /* -------------------------------------------------------------------------- */ template inline const Stored & ElementTypeMap:: operator()(const SupportType & type, const GhostType & ghost_type) const { auto it = this->getData(ghost_type).find(type); if (it == this->getData(ghost_type).end()) AKANTU_SILENT_EXCEPTION("No element of type " << ElementTypeMap::printType(type, ghost_type) << " in this ElementTypeMap<" << debug::demangle(typeid(Stored).name()) << "> class"); return it->second; } /* -------------------------------------------------------------------------- */ template inline Stored & ElementTypeMap:: operator()(const SupportType & type, const GhostType & ghost_type) { return this->getData(ghost_type)[type]; } /* -------------------------------------------------------------------------- */ template template inline Stored & ElementTypeMap:: operator()(U && insertee, const SupportType & type, const GhostType & ghost_type) { auto it = this->getData(ghost_type).find(type); if (it != this->getData(ghost_type).end()) { AKANTU_SILENT_EXCEPTION("Element of type " << ElementTypeMap::printType(type, ghost_type) << " already in this ElementTypeMap<" << debug::demangle(typeid(Stored).name()) << "> class"); } else { auto & data = this->getData(ghost_type); const auto & res = data.insert(std::make_pair(type, std::forward(insertee))); it = res.first; } return it->second; } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::DataMap & ElementTypeMap::getData(GhostType ghost_type) { if (ghost_type == _not_ghost) return data; return ghost_data; } /* -------------------------------------------------------------------------- */ template inline const typename ElementTypeMap::DataMap & ElementTypeMap::getData(GhostType ghost_type) const { if (ghost_type == _not_ghost) return data; return ghost_data; } /* -------------------------------------------------------------------------- */ /// Works only if stored is a pointer to a class with a printself method template void ElementTypeMap::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "ElementTypeMap<" << debug::demangle(typeid(Stored).name()) << "> [" << std::endl; for (auto gt : ghost_types) { const DataMap & data = getData(gt); for (auto & pair : data) { stream << space << space << ElementTypeMap::printType(pair.first, gt) << std::endl; } } stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ template ElementTypeMap::ElementTypeMap() = default; /* -------------------------------------------------------------------------- */ template ElementTypeMap::~ElementTypeMap() = default; /* -------------------------------------------------------------------------- */ /* ElementTypeMapArray */ /* -------------------------------------------------------------------------- */ template void ElementTypeMapArray::copy( const ElementTypeMapArray & other) { for (auto && ghost_type : ghost_types) { for (auto && type : this->elementTypes(_all_dimensions, ghost_types, _ek_not_defined)) { const auto & array_to_copy = other(type, ghost_type); auto & array = this->alloc(0, array_to_copy.getNbComponent(), type, ghost_type); array.copy(array_to_copy); } } } /* -------------------------------------------------------------------------- */ template inline Array & ElementTypeMapArray::alloc( UInt size, UInt nb_component, const SupportType & type, const GhostType & ghost_type, const T & default_value) { std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; Array * tmp; auto it = this->getData(ghost_type).find(type); if (it == this->getData(ghost_type).end()) { std::stringstream sstr; sstr << this->id << ":" << type << ghost_id; tmp = &(Memory::alloc(sstr.str(), size, nb_component, default_value)); std::stringstream sstrg; sstrg << ghost_type; // tmp->setTag(sstrg.str()); this->getData(ghost_type)[type] = tmp; } else { AKANTU_DEBUG_INFO( "The vector " << this->id << this->printType(type, ghost_type) << " already exists, it is resized instead of allocated."); tmp = it->second; it->second->resize(size); } return *tmp; } /* -------------------------------------------------------------------------- */ template inline void ElementTypeMapArray::alloc(UInt size, UInt nb_component, const SupportType & type, const T & default_value) { this->alloc(size, nb_component, type, _not_ghost, default_value); this->alloc(size, nb_component, type, _ghost, default_value); } /* -------------------------------------------------------------------------- */ template inline void ElementTypeMapArray::free() { AKANTU_DEBUG_IN(); for (auto gt : ghost_types) { auto & data = this->getData(gt); for (auto & pair : data) { dealloc(pair.second->getID()); } data.clear(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void ElementTypeMapArray::clear() { for (auto gt : ghost_types) { auto & data = this->getData(gt); for (auto & vect : data) { vect.second->clear(); } } } /* -------------------------------------------------------------------------- */ template template inline void ElementTypeMapArray::set(const ST & value) { for (auto gt : ghost_types) { auto & data = this->getData(gt); for (auto & vect : data) { vect.second->set(value); } } } /* -------------------------------------------------------------------------- */ template inline const Array & ElementTypeMapArray:: operator()(const SupportType & type, const GhostType & ghost_type) const { auto it = this->getData(ghost_type).find(type); if (it == this->getData(ghost_type).end()) AKANTU_SILENT_EXCEPTION("No element of type " << ElementTypeMapArray::printType(type, ghost_type) << " in this const ElementTypeMapArray<" << debug::demangle(typeid(T).name()) << "> class(\"" << this->id << "\")"); return *(it->second); } /* -------------------------------------------------------------------------- */ template inline Array & ElementTypeMapArray:: operator()(const SupportType & type, const GhostType & ghost_type) { auto it = this->getData(ghost_type).find(type); if (it == this->getData(ghost_type).end()) AKANTU_SILENT_EXCEPTION("No element of type " << ElementTypeMapArray::printType(type, ghost_type) << " in this ElementTypeMapArray<" << debug::demangle(typeid(T).name()) << "> class (\"" << this->id << "\")"); return *(it->second); } /* -------------------------------------------------------------------------- */ template inline void ElementTypeMapArray::setArray(const SupportType & type, const GhostType & ghost_type, const Array & vect) { auto it = this->getData(ghost_type).find(type); if (AKANTU_DEBUG_TEST(dblWarning) && it != this->getData(ghost_type).end() && it->second != &vect) { AKANTU_DEBUG_WARNING( "The Array " << this->printType(type, ghost_type) << " is already registred, this call can lead to a memory leak."); } this->getData(ghost_type)[type] = &(const_cast &>(vect)); } /* -------------------------------------------------------------------------- */ template inline void ElementTypeMapArray::onElementsRemoved( const ElementTypeMapArray & new_numbering) { for (auto gt : ghost_types) { for (auto & type : new_numbering.elementTypes(_all_dimensions, gt, _ek_not_defined)) { auto support_type = convertType(type); if (this->exists(support_type, gt)) { const auto & renumbering = new_numbering(type, gt); if (renumbering.size() == 0) continue; auto & vect = this->operator()(support_type, gt); auto nb_component = vect.getNbComponent(); Array tmp(renumbering.size(), nb_component); UInt new_size = 0; for (UInt i = 0; i < vect.size(); ++i) { UInt new_i = renumbering(i); if (new_i != UInt(-1)) { memcpy(tmp.storage() + new_i * nb_component, vect.storage() + i * nb_component, nb_component * sizeof(T)); ++new_size; } } tmp.resize(new_size); vect.copy(tmp); } } } } /* -------------------------------------------------------------------------- */ template void ElementTypeMapArray::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "ElementTypeMapArray<" << debug::demangle(typeid(T).name()) << "> [" << std::endl; for (UInt g = _not_ghost; g <= _ghost; ++g) { auto gt = (GhostType)g; const DataMap & data = this->getData(gt); typename DataMap::const_iterator it; for (it = data.begin(); it != data.end(); ++it) { stream << space << space << ElementTypeMapArray::printType(it->first, gt) << " [" << std::endl; it->second->printself(stream, indent + 3); stream << space << space << " ]" << std::endl; } } stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ /* SupportType Iterator */ /* -------------------------------------------------------------------------- */ template ElementTypeMap::type_iterator::type_iterator( DataMapIterator & list_begin, DataMapIterator & list_end, UInt dim, ElementKind ek) : list_begin(list_begin), list_end(list_end), dim(dim), kind(ek) {} /* -------------------------------------------------------------------------- */ template ElementTypeMap::type_iterator::type_iterator( const type_iterator & it) : list_begin(it.list_begin), list_end(it.list_end), dim(it.dim), kind(it.kind) {} /* -------------------------------------------------------------------------- */ template typename ElementTypeMap::type_iterator & ElementTypeMap::type_iterator:: operator=(const type_iterator & it) { if (this != &it) { list_begin = it.list_begin; list_end = it.list_end; dim = it.dim; kind = it.kind; } return *this; } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::type_iterator::reference ElementTypeMap::type_iterator::operator*() { return list_begin->first; } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::type_iterator::reference ElementTypeMap::type_iterator::operator*() const { return list_begin->first; } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::type_iterator & ElementTypeMap::type_iterator::operator++() { ++list_begin; while ((list_begin != list_end) && (((dim != _all_dimensions) && (dim != Mesh::getSpatialDimension(list_begin->first))) || ((kind != _ek_not_defined) && (kind != Mesh::getKind(list_begin->first))))) ++list_begin; return *this; } /* -------------------------------------------------------------------------- */ template typename ElementTypeMap::type_iterator ElementTypeMap::type_iterator::operator++(int) { type_iterator tmp(*this); operator++(); return tmp; } /* -------------------------------------------------------------------------- */ template inline bool ElementTypeMap::type_iterator:: operator==(const type_iterator & other) const { return this->list_begin == other.list_begin; } /* -------------------------------------------------------------------------- */ template inline bool ElementTypeMap::type_iterator:: operator!=(const type_iterator & other) const { return this->list_begin != other.list_begin; } /* -------------------------------------------------------------------------- */ template typename ElementTypeMap::ElementTypesIteratorHelper ElementTypeMap::elementTypesImpl(UInt dim, GhostType ghost_type, ElementKind kind) const { return ElementTypesIteratorHelper(*this, dim, ghost_type, kind); } /* -------------------------------------------------------------------------- */ template template typename ElementTypeMap::ElementTypesIteratorHelper ElementTypeMap::elementTypesImpl( const use_named_args_t & unused, pack &&... _pack) const { return ElementTypesIteratorHelper(*this, unused, _pack...); } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::type_iterator ElementTypeMap::firstType(UInt dim, GhostType ghost_type, ElementKind kind) const { typename DataMap::const_iterator b, e; b = getData(ghost_type).begin(); e = getData(ghost_type).end(); // loop until the first valid type while ((b != e) && (((dim != _all_dimensions) && (dim != Mesh::getSpatialDimension(b->first))) || ((kind != _ek_not_defined) && (kind != Mesh::getKind(b->first))))) ++b; return typename ElementTypeMap::type_iterator(b, e, dim, kind); } /* -------------------------------------------------------------------------- */ template inline typename ElementTypeMap::type_iterator ElementTypeMap::lastType(UInt dim, GhostType ghost_type, ElementKind kind) const { typename DataMap::const_iterator e; e = getData(ghost_type).end(); return typename ElementTypeMap::type_iterator(e, e, dim, kind); } /* -------------------------------------------------------------------------- */ /// standard output stream operator template inline std::ostream & operator<<(std::ostream & stream, const ElementTypeMap & _this) { _this.printself(stream); return stream; } /* -------------------------------------------------------------------------- */ -class ElementTypeMapArrayInializer { +class ElementTypeMapArrayInitializer { +protected: + using CompFunc = std::function; + public: - ElementTypeMapArrayInializer(UInt spatial_dimension = _all_dimensions, - UInt nb_component = 1, - const GhostType & ghost_type = _not_ghost, - const ElementKind & element_kind = _ek_regular) - : spatial_dimension(spatial_dimension), nb_component(nb_component), + ElementTypeMapArrayInitializer(const CompFunc & comp_func, + UInt spatial_dimension = _all_dimensions, + const GhostType & ghost_type = _not_ghost, + const ElementKind & element_kind = _ek_regular) + : comp_func(comp_func), spatial_dimension(spatial_dimension), ghost_type(ghost_type), element_kind(element_kind) {} const GhostType & ghostType() const { return ghost_type; } + virtual UInt nbComponent(const ElementType & type) const { + return comp_func(type, ghostType()); + } + protected: + CompFunc comp_func; UInt spatial_dimension; - UInt nb_component; GhostType ghost_type; ElementKind element_kind; }; /* -------------------------------------------------------------------------- */ -class MeshElementTypeMapArrayInializer : public ElementTypeMapArrayInializer { +class MeshElementTypeMapArrayInitializer + : public ElementTypeMapArrayInitializer { + using CompFunc = ElementTypeMapArrayInitializer::CompFunc; + public: - MeshElementTypeMapArrayInializer( + MeshElementTypeMapArrayInitializer( const Mesh & mesh, UInt nb_component = 1, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_regular, bool with_nb_element = false, bool with_nb_nodes_per_element = false) - : ElementTypeMapArrayInializer(spatial_dimension, nb_component, - ghost_type, element_kind), + : MeshElementTypeMapArrayInitializer( + mesh, + [nb_component](const ElementType &, const GhostType &) -> UInt { + return nb_component; + }, + spatial_dimension, ghost_type, element_kind, with_nb_element, + with_nb_nodes_per_element) {} + + MeshElementTypeMapArrayInitializer( + const Mesh & mesh, const CompFunc & comp_func, + UInt spatial_dimension = _all_dimensions, + const GhostType & ghost_type = _not_ghost, + const ElementKind & element_kind = _ek_regular, + bool with_nb_element = false, bool with_nb_nodes_per_element = false) + : ElementTypeMapArrayInitializer(comp_func, spatial_dimension, ghost_type, + element_kind), mesh(mesh), with_nb_element(with_nb_element), with_nb_nodes_per_element(with_nb_nodes_per_element) {} decltype(auto) elementTypes() const { return mesh.elementTypes(this->spatial_dimension, this->ghost_type, this->element_kind); } virtual UInt size(const ElementType & type) const { if (with_nb_element) return mesh.getNbElement(type, this->ghost_type); return 0; } - UInt nbComponent(const ElementType & type) const { + UInt nbComponent(const ElementType & type) const override { + UInt res = ElementTypeMapArrayInitializer::nbComponent(type); if (with_nb_nodes_per_element) - return (this->nb_component * mesh.getNbNodesPerElement(type)); + return (res * mesh.getNbNodesPerElement(type)); - return this->nb_component; + return res; } protected: const Mesh & mesh; bool with_nb_element; bool with_nb_nodes_per_element; }; /* -------------------------------------------------------------------------- */ -class FEEngineElementTypeMapArrayInializer - : public MeshElementTypeMapArrayInializer { +class FEEngineElementTypeMapArrayInitializer + : public MeshElementTypeMapArrayInitializer { public: - FEEngineElementTypeMapArrayInializer( + FEEngineElementTypeMapArrayInitializer( const FEEngine & fe_engine, UInt nb_component = 1, UInt spatial_dimension = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & element_kind = _ek_regular); + FEEngineElementTypeMapArrayInitializer( + const FEEngine & fe_engine, + const ElementTypeMapArrayInitializer::CompFunc & nb_component, + UInt spatial_dimension = _all_dimensions, + const GhostType & ghost_type = _not_ghost, + const ElementKind & element_kind = _ek_regular); + UInt size(const ElementType & type) const override; using ElementTypesIteratorHelper = ElementTypeMapArray::ElementTypesIteratorHelper; ElementTypesIteratorHelper elementTypes() const; protected: const FEEngine & fe_engine; }; /* -------------------------------------------------------------------------- */ template -template +template void ElementTypeMapArray::initialize(const Func & f, const T & default_value, - bool do_not_default, - CompFunc && comp_func) { + bool do_not_default) { auto ghost_type = f.ghostType(); for (auto & type : f.elementTypes()) { if (not this->exists(type, ghost_type)) if (do_not_default) { - auto & array = - this->alloc(0, comp_func(type, ghost_type), type, ghost_type); + auto & array = this->alloc(0, f.nbComponent(type), type, ghost_type); array.resize(f.size(type)); } else { - this->alloc(f.size(type), comp_func(type, ghost_type), type, ghost_type, + this->alloc(f.size(type), f.nbComponent(type), type, ghost_type, default_value); } else { auto & array = this->operator()(type, ghost_type); if (not do_not_default) array.resize(f.size(type), default_value); else array.resize(f.size(type)); } } } /* -------------------------------------------------------------------------- */ +/** + * All parameters are named optionals + * \param _nb_component a functor giving the number of components per + * (ElementType, GhostType) pair or a scalar giving a unique number of + * components + * regardless of type + * \param _spatial_dimension a filter for the elements of a specific dimension + * \param _element_kind filter with element kind (_ek_regular, _ek_structural, + * ...) + * \param _with_nb_element allocate the arrays with the number of elements for + * each + * type in the mesh + * \param _with_nb_nodes_per_element multiply the number of components by the + * number of nodes per element + * \param _default_value default inital value + * \param _do_not_default do not initialize the allocated arrays + * \param _ghost_type filter a type of ghost + */ template template void ElementTypeMapArray::initialize(const Mesh & mesh, pack &&... _pack) { GhostType requested_ghost_type = OPTIONAL_NAMED_ARG(ghost_type, _casper); bool all_ghost_types = requested_ghost_type == _casper; for (auto ghost_type : ghost_types) { if ((not(ghost_type == requested_ghost_type)) and (not all_ghost_types)) continue; - auto functor = MeshElementTypeMapArrayInializer( - mesh, OPTIONAL_NAMED_ARG(nb_component, 1), + // This thing should have a UInt or functor type + auto && nb_components = OPTIONAL_NAMED_ARG(nb_component, 1); + auto functor = MeshElementTypeMapArrayInitializer( + mesh, std::forward(nb_components), OPTIONAL_NAMED_ARG(spatial_dimension, mesh.getSpatialDimension()), ghost_type, OPTIONAL_NAMED_ARG(element_kind, _ek_regular), OPTIONAL_NAMED_ARG(with_nb_element, false), OPTIONAL_NAMED_ARG(with_nb_nodes_per_element, false)); - std::function - nb_component_functor = - [&](const ElementType & type, const GhostType & - /*ghost_type*/) -> UInt { return functor.nbComponent(type); }; - - this->initialize( - functor, OPTIONAL_NAMED_ARG(default_value, T()), - OPTIONAL_NAMED_ARG(do_not_default, false), - OPTIONAL_NAMED_ARG(nb_component_functor, nb_component_functor)); + this->initialize(functor, OPTIONAL_NAMED_ARG(default_value, T()), + OPTIONAL_NAMED_ARG(do_not_default, false)); } } /* -------------------------------------------------------------------------- */ +/** + * All parameters are named optionals + * \param _nb_component a functor giving the number of components per + * (ElementType, GhostType) pair or a scalar giving a unique number of + * components + * regardless of type + * \param _spatial_dimension a filter for the elements of a specific dimension + * \param _element_kind filter with element kind (_ek_regular, _ek_structural, + * ...) + * \param _default_value default inital value + * \param _do_not_default do not initialize the allocated arrays + * \param _ghost_type filter a specific ghost type + * \param _all_ghost_types get all ghost types + */ template template void ElementTypeMapArray::initialize(const FEEngine & fe_engine, pack &&... _pack) { bool all_ghost_types = OPTIONAL_NAMED_ARG(all_ghost_types, true); GhostType requested_ghost_type = OPTIONAL_NAMED_ARG(ghost_type, _not_ghost); for (auto ghost_type : ghost_types) { if ((not(ghost_type == requested_ghost_type)) and (not all_ghost_types)) continue; - auto functor = FEEngineElementTypeMapArrayInializer( - fe_engine, OPTIONAL_NAMED_ARG(nb_component, 1), + auto && nb_components = OPTIONAL_NAMED_ARG(nb_component, 1); + auto functor = FEEngineElementTypeMapArrayInitializer( + fe_engine, std::forward(nb_components), OPTIONAL_NAMED_ARG(spatial_dimension, UInt(-2)), ghost_type, OPTIONAL_NAMED_ARG(element_kind, _ek_regular)); - std::function - nb_component_functor = - [&](const ElementType & type, const GhostType & - /*ghost_type*/) -> UInt { return functor.nbComponent(type); }; - - this->initialize( - functor, OPTIONAL_NAMED_ARG(default_value, T()), - OPTIONAL_NAMED_ARG(do_not_default, false), - OPTIONAL_NAMED_ARG(nb_component_functor, nb_component_functor)); + this->initialize(functor, OPTIONAL_NAMED_ARG(default_value, T()), + OPTIONAL_NAMED_ARG(do_not_default, false)); } } /* -------------------------------------------------------------------------- */ template inline T & ElementTypeMapArray:: operator()(const Element & element, UInt component) { return this->operator()(element.type, element.ghost_type)(element.element, component); } /* -------------------------------------------------------------------------- */ template inline const T & ElementTypeMapArray:: operator()(const Element & element, UInt component) const { return this->operator()(element.type, element.ghost_type)(element.element, component); } } // namespace akantu #endif /* __AKANTU_ELEMENT_TYPE_MAP_TMPL_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 09caa4037..4dc96d2e7 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,818 +1,803 @@ /** * @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), 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(); 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); + foreground->initialize(emodel.getMesh(), _nb_component = 1, + _ghost_type = gt); + background->initialize(emodel.getMesh(), _nb_component = 1, + _ghost_type = gt); // 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 { + class BackgroundShapeDInitializer : public ElementTypeMapArrayInitializer { public: - BackgroundShapeDInitializer(UInt spatial_dimension, - FEEngine & engine, - const ElementType & foreground_type, - ElementTypeMapArray & filter, + BackgroundShapeDInitializer(UInt spatial_dimension, FEEngine & engine, + const ElementType & foreground_type, + const ElementTypeMapArray & filter, const GhostType & ghost_type) - : ElementTypeMapArrayInializer(spatial_dimension, 0, ghost_type, - _ek_regular) { + : ElementTypeMapArrayInitializer( + [](const ElementType & bgtype, const GhostType &) { + return ShapeFunctions::getShapeDerivativesSize(bgtype); + }, + spatial_dimension, 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/src/model/structural_mechanics/structural_elements/structural_element_kirchhoff_shell.hh b/src/model/structural_mechanics/structural_elements/structural_element_kirchhoff_shell.hh index 6e9806e72..68a85eac0 100644 --- a/src/model/structural_mechanics/structural_elements/structural_element_kirchhoff_shell.hh +++ b/src/model/structural_mechanics/structural_elements/structural_element_kirchhoff_shell.hh @@ -1,75 +1,76 @@ /** * @file structural_element_bernoulli_kirchhoff_shell.hh * * @author Fabian Barras * @author Sébastien Hartmann * @author Nicolas Richart * @author Damien Spielmann * * @date creation Tue Sep 19 2017 * * @brief Specific functions for bernoulli kirchhoff shell * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_STRUCTURAL_ELEMENT_BERNOULLI_KIRCHHOFF_SHELL_HH__ #define __AKANTU_STRUCTURAL_ELEMENT_BERNOULLI_KIRCHHOFF_SHELL_HH__ #include "structural_mechanics_model.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template <> inline void StructuralMechanicsModel::assembleMass<_discrete_kirchhoff_triangle_18>() { AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ template <> void StructuralMechanicsModel::computeTangentModuli< _discrete_kirchhoff_triangle_18>(Array & tangent_moduli) { auto tangent_size = ElementClass<_discrete_kirchhoff_triangle_18>::getNbStressComponents(); auto nb_quad = getFEEngine().getNbIntegrationPoints(_discrete_kirchhoff_triangle_18); auto H_it = tangent_moduli.begin(tangent_size, tangent_size); for (UInt mat : element_material(_discrete_kirchhoff_triangle_18, _not_ghost)) { auto & m = materials[mat]; for (UInt q = 0; q < nb_quad; ++q, ++H_it) { auto & H = *H_it; + H.clear(); Matrix D = {{1, m.nu, 0}, {m.nu, 1, 0}, {0, 0, (1 - m.nu) / 2}}; D *= m.E / (1 - m.nu * m.nu); - H.block(0, 0, 3, 3) = D; // in plane membrane behavior - H.block(3, 3, 3, 3) = D * Math::pow<3>(m.t) / 12.; // bending behavior + H.block(D, 0, 0); // in plane membrane behavior + H.block(D * Math::pow<3>(m.t) / 12., 3, 3); // bending behavior } } } } // akantu #endif /* __AKANTU_STRUCTURAL_ELEMENT_BERNOULLI_DISCRETE_KIRCHHOFF_TRIANGLE_18_HH__ */ diff --git a/src/model/structural_mechanics/structural_mechanics_model.cc b/src/model/structural_mechanics/structural_mechanics_model.cc index dc3e5e8d2..8ed981deb 100644 --- a/src/model/structural_mechanics/structural_mechanics_model.cc +++ b/src/model/structural_mechanics/structural_mechanics_model.cc @@ -1,442 +1,490 @@ /** * @file structural_mechanics_model.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 Model implementation for StucturalMechanics elements * * @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" #include "dof_manager.hh" #include "integrator_gauss.hh" #include "mesh.hh" #include "shape_structural.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER #include "dumpable_inline_impl.hh" #include "dumper_elemental_field.hh" #include "dumper_iohelper_paraview.hh" #include "group_manager_inline_impl.cc" #endif /* -------------------------------------------------------------------------- */ #include "structural_mechanics_model_inline_impl.cc" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ StructuralMechanicsModel::StructuralMechanicsModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : Model(mesh, ModelType::_structural_mechanics_model, dim, id, memory_id), time_step(NAN), 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_DEBUG_TO_IMPLEMENT(); } #ifdef AKANTU_USE_IOHELPER this->mesh.registerDumper("paraview_all", id, true); #endif this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_structural); this->initDOFManager(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ StructuralMechanicsModel::~StructuralMechanicsModel() = default; /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initFullImpl(const ModelOptions & options) { // <<<< This is the SolidMechanicsModel implementation for future ref >>>> // material_index.initialize(mesh, _element_kind = _ek_not_defined, // _default_value = UInt(-1), _with_nb_element = // true); // material_local_numbering.initialize(mesh, _element_kind = _ek_not_defined, // _with_nb_element = true); // Model::initFullImpl(options); // // initialize pbc // if (this->pbc_pair.size() != 0) // this->initPBC(); // // initialize the materials // if (this->parser.getLastParsedFile() != "") { // this->instantiateMaterials(); // } // this->initMaterials(); // this->initBC(*this, *displacement, *displacement_increment, // *external_force); // <<<< END >>>> 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, + _all_ghost_types = true, + _nb_component = [&stress_components]( + const ElementType & type, const GhostType &) -> 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) { // this->time_step = time_step; // #if defined(AKANTU_USE_IOHELPER) // this->mesh.getDumper().setTimeStep(time_step); // #endif // } /* -------------------------------------------------------------------------- */ /* Initialisation */ /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::initSolver( TimeStepSolverType time_step_solver_type, NonLinearSolverType) { AKANTU_DEBUG_IN(); this->allocNodalField(displacement_rotation, nb_degree_of_freedom, "displacement"); - this->allocNodalField(external_force_momentum, nb_degree_of_freedom, + this->allocNodalField(external_force, nb_degree_of_freedom, "external_force"); - this->allocNodalField(internal_force_momentum, nb_degree_of_freedom, + 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 (!dof_manager.hasDOFs("displacement")) { dof_manager.registerDOFs("displacement", *displacement_rotation, _dst_nodal); dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs); } if (time_step_solver_type == _tsst_dynamic || time_step_solver_type == _tsst_dynamic_lumped) { this->allocNodalField(velocity, spatial_dimension, "velocity"); this->allocNodalField(acceleration, spatial_dimension, "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() { for (auto && type : mesh.elementTypes(_element_kind = _ek_structural)) { // computeRotationMatrix(type); element_material.alloc(mesh.getNbElement(type), 1, type); } getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); getDOFManager().getMatrix("K").clear(); for (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 } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::computeStresses() { AKANTU_DEBUG_IN(); for (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(); } -/* -------------------------------------------------------------------------- */ -void StructuralMechanicsModel::updateResidual() { - AKANTU_DEBUG_IN(); - - auto & K = getDOFManager().getMatrix("K"); - - internal_force_momentum->clear(); - K.matVecMul(*displacement_rotation, *internal_force_momentum); - *internal_force_momentum *= -1.; - - getDOFManager().assembleToResidual("displacement", *external_force_momentum, - 1.); - - getDOFManager().assembleToResidual("displacement", *internal_force_momentum, - 1.); - - AKANTU_DEBUG_OUT(); -} - /* -------------------------------------------------------------------------- */ void StructuralMechanicsModel::computeRotationMatrix(const ElementType & type) { Mesh & mesh = getFEEngine().getMesh(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_element = mesh.getNbElement(type); if (!rotation_matrix.exists(type)) { rotation_matrix.alloc(nb_element, nb_degree_of_freedom * nb_nodes_per_element * nb_degree_of_freedom * nb_nodes_per_element, type); } else { rotation_matrix(type).resize(nb_element); } rotation_matrix(type).clear(); Array rotations(nb_element, nb_degree_of_freedom * nb_degree_of_freedom); rotations.clear(); #define COMPUTE_ROTATION_MATRIX(type) computeRotationMatrix(rotations); AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(COMPUTE_ROTATION_MATRIX); #undef COMPUTE_ROTATION_MATRIX auto R_it = rotations.begin(nb_degree_of_freedom, nb_degree_of_freedom); auto T_it = rotation_matrix(type).begin(nb_degree_of_freedom * nb_nodes_per_element, nb_degree_of_freedom * nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++R_it, ++T_it) { auto & T = *T_it; auto & R = *R_it; for (UInt k = 0; k < nb_nodes_per_element; ++k) { for (UInt i = 0; i < nb_degree_of_freedom; ++i) for (UInt j = 0; j < nb_degree_of_freedom; ++j) T(k * nb_degree_of_freedom + i, k * nb_degree_of_freedom + j) = R(i, j); } } } /* -------------------------------------------------------------------------- */ dumper::Field * 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; dumper::Field * field = NULL; field = mesh.createNodalField(uint_nodal_fields[field_name], group_name); return field; } /* -------------------------------------------------------------------------- */ dumper::Field * 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; } if (field_name == "displacement") { return mesh.createStridedNodalField(displacement_rotation, group_name, n, 0, padding_flag); } if (field_name == "rotation") { return mesh.createStridedNodalField(displacement_rotation, group_name, nb_degree_of_freedom - n, n, padding_flag); } if (field_name == "force") { - return mesh.createStridedNodalField(external_force_momentum, group_name, n, + return mesh.createStridedNodalField(external_force, group_name, n, 0, padding_flag); } if (field_name == "momentum") { - return mesh.createStridedNodalField(external_force_momentum, group_name, + return mesh.createStridedNodalField(external_force, group_name, nb_degree_of_freedom - n, n, padding_flag); } if (field_name == "internal_force") { - return mesh.createStridedNodalField(internal_force_momentum, group_name, n, + return mesh.createStridedNodalField(internal_force, group_name, n, 0, padding_flag); } if (field_name == "internal_momentum") { - return mesh.createStridedNodalField(internal_force_momentum, group_name, + return mesh.createStridedNodalField(internal_force, group_name, nb_degree_of_freedom - n, n, padding_flag); } return nullptr; } /* -------------------------------------------------------------------------- */ dumper::Field * StructuralMechanicsModel::createElementalField( const std::string & field_name, const std::string & group_name, bool, const ElementKind & kind, const std::string &) { dumper::Field * field = NULL; if (field_name == "element_index_by_material") field = mesh.createElementalField( field_name, group_name, this->spatial_dimension, kind); 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(); } /// callback to assemble a lumped Matrix void StructuralMechanicsModel::assembleLumpedMatrix(const ID & /*id*/) {} /// callback to assemble the residual StructuralMechanicsModel::(rhs) void StructuralMechanicsModel::assembleResidual() { - this->getDOFManager().assembleToResidual("displacement", - *this->external_force_momentum, 1); + AKANTU_DEBUG_IN(); + + internal_force->clear(); + computeStresses(); + assembleInternalForce(); + + getDOFManager().assembleToResidual("displacement", *external_force, + 1.); + + getDOFManager().assembleToResidual("displacement", *internal_force, + -1); + + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Virtual methods from MeshEventHandler */ /* -------------------------------------------------------------------------- */ /// function to implement to react on akantu::NewNodesEvent void StructuralMechanicsModel::onNodesAdded(const Array & /*nodes_list*/, const NewNodesEvent & /*event*/) {} /// function to implement to react on akantu::RemovedNodesEvent void StructuralMechanicsModel::onNodesRemoved( const Array & /*nodes_list*/, const Array & /*new_numbering*/, const RemovedNodesEvent & /*event*/) {} /// function to implement to react on akantu::NewElementsEvent void StructuralMechanicsModel::onElementsAdded( const Array & /*elements_list*/, const NewElementsEvent & /*event*/) {} /// function to implement to react on akantu::RemovedElementsEvent void StructuralMechanicsModel::onElementsRemoved( const Array & /*elements_list*/, const ElementTypeMapArray & /*new_numbering*/, const RemovedElementsEvent & /*event*/) {} /// function to implement to react on akantu::ChangedElementsEvent void StructuralMechanicsModel::onElementsChanged( const Array & /*old_elements_list*/, const Array & /*new_elements_list*/, const ElementTypeMapArray & /*new_numbering*/, const ChangedElementsEvent & /*event*/) {} /* -------------------------------------------------------------------------- */ /* 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", _tsst_static); } case _implicit_dynamic: { return std::make_tuple("implicit", _tsst_dynamic); } default: return std::make_tuple("unknown", _tsst_not_defined); } } /* ------------------------------------------------------------------------ */ ModelSolverOptions StructuralMechanicsModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case _tsst_static: { options.non_linear_solver_type = _nls_linear; options.integration_scheme_type["displacement"] = _ist_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } default: AKANTU_EXCEPTION(type << " is not a valid time step solver type"); } return options; } + +/* -------------------------------------------------------------------------- */ +void StructuralMechanicsModel::assembleInternalForce() { + for (auto type : mesh.elementTypes(_spatial_dimension = _all_dimensions, + _element_kind = _ek_structural)) { + assembleInternalForce(type, _not_ghost); + // assembleInternalForce(type, _ghost); + } +} + +/* -------------------------------------------------------------------------- */ +void StructuralMechanicsModel::assembleInternalForce(const 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 fcb036768..b06e1ea41 100644 --- a/src/model/structural_mechanics/structural_mechanics_model.hh +++ b/src/model/structural_mechanics/structural_mechanics_model.hh @@ -1,327 +1,333 @@ /** * @file structural_mechanics_model.hh * * @author Fabian Barras * @author Sébastien Hartmann * @author Nicolas Richart * @author Damien Spielmann * * @date creation: Fri Jul 15 2011 * @date last modification: Thu Jan 21 2016 * * @brief Particular implementation of the structural elements in the * 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 "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 spatial_dimension = _all_dimensions, const ID & id = "structural_mechanics_model", const MemoryID & memory_id = 0); virtual ~StructuralMechanicsModel(); /// 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 &) override; /// callback to assemble a Matrix void assembleMatrix(const ID &) override; /// callback to assemble a lumped Matrix void assembleLumpedMatrix(const ID &) override; /// callback to assemble the residual (rhs) void assembleResidual() override; /* ------------------------------------------------------------------------ */ /* Virtual methods from MeshEventHandler */ /* ------------------------------------------------------------------------ */ /// function to implement to react on akantu::NewNodesEvent void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) override; /// function to implement to react on akantu::RemovedNodesEvent void onNodesRemoved(const Array & nodes_list, const Array & new_numbering, const RemovedNodesEvent & event) override; /// function to implement to react on akantu::NewElementsEvent void onElementsAdded(const Array & elements_list, const NewElementsEvent & event) override; /// function to implement to react on akantu::RemovedElementsEvent void onElementsRemoved(const Array & elements_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) override; /// function to implement to react on akantu::ChangedElementsEvent void onElementsChanged(const Array & old_elements_list, const Array & new_elements_list, const ElementTypeMapArray & new_numbering, const ChangedElementsEvent & event) override; /* ------------------------------------------------------------------------ */ /* 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; + UInt getNbDegreeOfFreedom(const ElementType & type) const; + /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ void initSolver(TimeStepSolverType, NonLinearSolverType) 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(const ElementType & type, GhostType gt); + /// assemble the stiffness matrix void assembleStiffnessMatrix(); /// assemble the mass matrix for consistent mass resolutions void assembleMass(); - /// update the residual vector - void updateResidual(); - + /// TODO remove void computeRotationMatrix(const ElementType & type); protected: /// compute Rotation Matrices template void computeRotationMatrix(__attribute__((unused)) Array & rotations) {} /* ------------------------------------------------------------------------ */ /* Mass (structural_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMass(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 assembleMass(); template void computeStressOnQuad(); template void computeTangentModuli(Array & tangent_moduli); /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: virtual dumper::Field * createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag); virtual dumper::Field * createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag); virtual dumper::Field * createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, const ElementKind & kind, const std::string & fe_engine_id = ""); /* ------------------------------------------------------------------------ */ /* 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_momentum, Array &); + AKANTU_GET_MACRO(ExternalForce, *external_force, Array &); /// get the StructuralMechanicsModel::internal_force vector (boundary forces) - AKANTU_GET_MACRO(InternalForce, *internal_force_momentum, Array &); + 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); void addMaterial(StructuralMaterial & material) { materials.push_back(material); } const StructuralMaterial & getMaterial(const Element & element) const { return materials[element_material(element)]; } /* ------------------------------------------------------------------------ */ /* 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_momentum{nullptr}; + Array * internal_force{nullptr}; /// forces array - Array * external_force_momentum{nullptr}; + 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; /* ------------------------------------------------------------------------ */ std::vector materials; }; } // namespace akantu #endif /* __AKANTU_STRUCTURAL_MECHANICS_MODEL_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 ec27c7ede..7ef7ed0cc 100644 --- a/src/model/structural_mechanics/structural_mechanics_model_inline_impl.cc +++ b/src/model/structural_mechanics/structural_mechanics_model_inline_impl.cc @@ -1,361 +1,370 @@ /** * @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 { +/* -------------------------------------------------------------------------- */ +inline UInt +StructuralMechanicsModel::getNbDegreeOfFreedom(const ElementType & type) const { + UInt ndof = 0; +#define GET_(type) ndof = ElementClass::getNbDegreeOfFreedom() + AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_, _ek_structural); +#undef GET_ + + return ndof; +} + /* -------------------------------------------------------------------------- */ 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 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); + // No need to rotate because B is post-multiplied + for (UInt e = 0; e < nb_element; ++e, ++ug) { for (UInt q = 0; q < nb_quadrature_points; ++q, ++D_B, ++DBu) { - DBu->template mul(*D_B, ul); + DBu->template mul(*D_B, *ug); } } 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__ */ 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 b2d1070d9..f33fd4e01 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,99 +1,107 @@ /** * @file test_structural_mechanics_model_bernoulli_beam_2.cc * * @author Fabian Barras * @author Lucas Frérot * * @date creation: Fri Jul 15 2011 * @date last modification: Sun Oct 19 2014 * * @brief Computation of the analytical exemple 1.1 in the TGC vol 6 * * @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 "test_structural_mechanics_model_fixture.hh" +#include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ #include using namespace akantu; /* -------------------------------------------------------------------------- */ class TestStructDKT18 : public TestStructuralFixture> { 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 { 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().clear(); 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 - *disp = {40, 20, -800 , -20, 40, 0}; ++disp; - *disp = {50, 40, -1400, -40, 50, 0}; ++disp; - *disp = {10, 20, -200 , -20, 10, 0}; ++disp; + // *disp = {40, 20, -800 , -20, 40, 0}; ++disp; + // *disp = {50, 40, -1400, -40, 50, 0}; ++disp; + // *disp = {10, 20, -200 , -20, 10, 0}; ++disp; + *disp = {0, 0, -800 , -20, 40, 0}; ++disp; + *disp = {0, 0, -1400, -40, 50, 0}; ++disp; + *disp = {0, 0, -200 , -20, 10, 0}; ++disp; // clang-format on } void setNeumanns() override {} protected: StructuralMaterial mat; }; /* -------------------------------------------------------------------------- */ // Batoz Vol 2. patch test, ISBN 2-86601-259-3 TEST_F(TestStructDKT18, TestDisplacements) { model->solveStep(); + model->getDOFManager().getMatrix("K").saveMatrix("K.mtx"); + model->getDOFManager().getMatrix("J").saveMatrix("J.mtx"); 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; std::cout << center_node_disp << std::endl; + auto error = solution - center_node_disp; + Real tol = Math::getTolerance(); EXPECT_NEAR(error.norm(), 0., tol); }