diff --git a/src/fe_engine/shape_structural.hh b/src/fe_engine/shape_structural.hh index 057326226..58bfef36e 100644 --- a/src/fe_engine/shape_structural.hh +++ b/src/fe_engine/shape_structural.hh @@ -1,120 +1,136 @@ /** * @file shape_structural.hh * * @author Fabian Barras * @author Nicolas Richart * * @date creation: Tue Feb 15 2011 * @date last modification: Thu Oct 22 2015 * * @brief shape class for element with different set of shapes functions * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "shape_functions.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SHAPE_STRUCTURAL_HH__ #define __AKANTU_SHAPE_STRUCTURAL_HH__ namespace akantu { template class ShapeStructural : public ShapeFunctions { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: using ElementTypeMapMultiReal = ElementTypeMap **>; ShapeStructural(Mesh & mesh, const ID & id = "shape_structural", const MemoryID & memory_id = 0); ~ShapeStructural() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// compute shape functions on given integration points template void computeShapesOnIntegrationPoints( const Array &, const Matrix & integration_points, Array & shapes, const GhostType & ghost_type, const Array & filter_elements = empty_filter) const; /// initialization function for structural elements inline void initShapeFunctions(const Array & nodes, const Matrix & integration_points, const ElementType & type, const GhostType & ghost_type); /// precompute the rotation matrices for the elements dofs template void precomputeRotationMatrices(const Array & nodes, const GhostType & ghost_type); /// pre compute all shapes on the element integration points from natural /// coordinates template void precomputeShapesOnIntegrationPoints(const Array & nodes, const GhostType & ghost_type); /// pre compute all shapes on the element integration points from natural /// coordinates template void precomputeShapeDerivativesOnIntegrationPoints(const Array & nodes, const GhostType & ghost_type); /// interpolate nodal values on the integration points template void interpolateOnIntegrationPoints( const Array & u, Array & uq, UInt nb_degree_of_freedom, const GhostType & ghost_type = _not_ghost, const Array & filter_elements = empty_filter) const; /// compute the gradient of u on the integration points template void gradientOnIntegrationPoints( const Array & u, Array & nablauq, UInt nb_degree_of_freedom, const GhostType & ghost_type = _not_ghost, const Array & filter_elements = empty_filter) const; + /// interpolate on physical point + template + void interpolate(const Vector & /*real_coords*/, UInt /*elem*/, + const Matrix & /*nodal_values*/, + Vector & /*interpolated*/, + const GhostType & /*ghost_type*/) const { + AKANTU_DEBUG_TO_IMPLEMENT(); + } + + template + void computeShapes(const Vector & /*real_coords*/, UInt /*elem*/, + Vector & /*shapes*/, + const GhostType & /*ghost_type*/) const { + AKANTU_DEBUG_TO_IMPLEMENT(); + } + /// multiply a field by shape functions template void fieldTimesShapes(__attribute__((unused)) const Array & field, - __attribute__((unused)) Array & fiedl_times_shapes, + __attribute__((unused)) Array & field_times_shapes, __attribute__((unused)) const GhostType & ghost_type) const { AKANTU_DEBUG_TO_IMPLEMENT(); } protected: ElementTypeMapArray rotation_matrices; }; } // namespace akantu #include "shape_structural_inline_impl.cc" #endif /* __AKANTU_SHAPE_STRUCTURAL_HH__ */ diff --git a/src/fe_engine/shape_structural_inline_impl.cc b/src/fe_engine/shape_structural_inline_impl.cc index 3c39cb280..345b0f519 100644 --- a/src/fe_engine/shape_structural_inline_impl.cc +++ b/src/fe_engine/shape_structural_inline_impl.cc @@ -1,353 +1,353 @@ /** * @file shape_structural_inline_impl.cc * * @author Fabian Barras * @author Nicolas Richart * * @date creation: Mon Dec 13 2010 * @date last modification: Thu Oct 15 2015 * * @brief ShapeStructural inline implementation * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "mesh_iterators.hh" #include "shape_structural.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SHAPE_STRUCTURAL_INLINE_IMPL_CC__ #define __AKANTU_SHAPE_STRUCTURAL_INLINE_IMPL_CC__ namespace akantu { template inline void ShapeStructural::initShapeFunctions( const Array & /* unused */, const Matrix & /* unused */, const ElementType & /* unused */, const GhostType & /* unused */) { AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ #define INIT_SHAPE_FUNCTIONS(type) \ setIntegrationPointsByType(integration_points, ghost_type); \ 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 { /// \TODO this code differs from ShapeLagrangeBase only in the size of N UInt nb_points = integration_points.cols(); UInt nb_element = mesh.getConnectivity(type, ghost_type).size(); shapes.resize(nb_element * nb_points); UInt ndof = ElementClass::getNbDegreeOfFreedom(); #if !defined(AKANTU_NDEBUG) UInt size_of_shapes = ElementClass::getShapeSize(); AKANTU_DEBUG_ASSERT(shapes.getNbComponent() == size_of_shapes, "The shapes array does not have the correct " << "number of component"); #endif auto shapes_it = shapes.begin_reinterpret( ElementClass::getNbNodesPerInterpolationElement(), ndof, nb_points, nb_element); auto shapes_begin = shapes_it; if (filter_elements != empty_filter) { nb_element = filter_elements.size(); } for (UInt elem = 0; elem < nb_element; ++elem) { if (filter_elements != empty_filter) shapes_it = shapes_begin + filter_elements(elem); Tensor3 & N = *shapes_it; ElementClass::computeShapes(integration_points, N); if (filter_elements == empty_filter) ++shapes_it; } } /* -------------------------------------------------------------------------- */ template template void ShapeStructural::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::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()); } } 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(); auto itp_type = FEEngine::getInterpolationType(type); if (not shapes.exists(itp_type, ghost_type)) { auto size_of_shapes = this->getShapeSize(type); this->shapes.alloc(0, size_of_shapes, itp_type, ghost_type); } auto & shapes_ = this->shapes(itp_type, ghost_type); shapes_.resize(nb_element * nb_points); auto shapes_it = shapes_.begin_reinterpret( nb_dof, nb_dof * nb_nodes_per_element, nb_points, nb_element); for (UInt elem = 0; elem < nb_element; ++elem, ++shapes_it) { auto & N = *shapes_it; ElementClass::computeShapes(natural_coords, N); } AKANTU_DEBUG_OUT(); } // namespace akantu /* -------------------------------------------------------------------------- */ template template void ShapeStructural::precomputeShapeDerivativesOnIntegrationPoints( const Array & nodes, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); const auto & natural_coords = integration_points(type, ghost_type); 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); auto R = RDOFs.block(0, 0, spatial_dimension, spatial_dimension); // Rotate to local basis auto x = (R * X).block(0, 0, natural_spatial_dimension, nb_nodes_per_element); Tensor3 dnds(B.size(0), B.size(1), B.size(2)); ElementClass::computeDNDS(natural_coords, dnds); Tensor3 J(x.rows(), natural_coords.rows(), natural_coords.cols()); ElementClass::computeJMat(dnds, x, J); ElementClass::computeShapeDerivatives(J, dnds, B); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void ShapeStructural::interpolateOnIntegrationPoints( const Array & in_u, Array & out_uq, UInt nb_dof, const GhostType & ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(out_uq.getNbComponent() == nb_dof, "The output array shape is not correct"); auto itp_type = FEEngine::getInterpolationType(type); const auto & shapes_ = shapes(itp_type, ghost_type); auto nb_element = mesh.getNbElement(type, ghost_type); auto nb_nodes_per_element = ElementClass::getNbNodesPerElement(); auto nb_quad_points_per_element = integration_points(type, ghost_type).cols(); Array u_el(0, nb_nodes_per_element * nb_dof); - FEEngine::extractNodalToElementField(mesh, in_u, u_el, ghost_type, - filter_elements); + 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_elements(nb_element, filter_elements, [&](auto && el) { auto & uq = *out_it; const auto & u = *u_it; auto N = Tensor3(shapes_it[el]); for (auto && q : arange(uq.size(2))) { auto uq_q = Matrix(uq(q)); auto u_q = Matrix(u(q)); auto N_q = Matrix(N(q)); uq_q.mul(N, u); } ++out_it; ++u_it; }); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template void ShapeStructural::gradientOnIntegrationPoints( const Array & in_u, Array & out_nablauq, UInt nb_dof, const GhostType & ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); auto itp_type = FEEngine::getInterpolationType(type); const auto & shapesd = shapes_derivatives(itp_type, ghost_type); auto nb_element = mesh.getNbElement(type, ghost_type); auto element_dimension = ElementClass::getSpatialDimension(); auto nb_quad_points_per_element = integration_points(type, ghost_type).cols(); auto nb_nodes_per_element = ElementClass::getNbNodesPerElement(); Array u_el(0, nb_nodes_per_element * nb_dof); FEEngine::extractNodalToElementField(mesh, in_u, u_el, ghost_type, filter_elements); auto nb_quad_points = nb_quad_points_per_element * u_el.size(); out_nablauq.resize(nb_quad_points); auto out_it = out_nablauq.begin_reinterpret( element_dimension, 1, nb_quad_points_per_element, u_el.size()); auto shapesd_it = shapesd.begin_reinterpret( element_dimension, nb_dof * nb_nodes_per_element, nb_quad_points_per_element, nb_element); auto u_it = u_el.begin_reinterpret(nb_dof * nb_nodes_per_element, 1, nb_quad_points_per_element, u_el.size()); for_each_elements(nb_element, filter_elements, [&](auto && el) { auto & nablau = *out_it; const auto & u = *u_it; auto B = Tensor3(shapesd_it[el]); for (auto && q : arange(nablau.size(2))) { auto nablau_q = Matrix(uq(q)); auto u_q = Matrix(u(q)); auto B_q = Matrix(N(q)); nablau_q.mul(B, u); } ++out_it; ++u_it; }); AKANTU_DEBUG_OUT(); } } // namespace akantu #endif /* __AKANTU_SHAPE_STRUCTURAL_INLINE_IMPL_CC__ */ diff --git a/test/test_fe_engine/CMakeLists.txt b/test/test_fe_engine/CMakeLists.txt index 318d1a315..653f0b09d 100644 --- a/test/test_fe_engine/CMakeLists.txt +++ b/test/test_fe_engine/CMakeLists.txt @@ -1,113 +1,113 @@ #=============================================================================== # @file CMakeLists.txt # # @author Guillaume Anciaux # # @date creation: Fri Sep 03 2010 # @date last modification: Mon Dec 07 2015 # # @brief configuration for FEM tests # # @section LICENSE # # Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de # Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des # Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== #=============================================================================== function(register_fem_test operation type) set(_target test_${operation}${type}) register_test(${_target} SOURCES test_${operation}.cc FILES_TO_COPY ${type}.msh COMPILE_OPTIONS TYPE=${type} PACKAGE core ) endfunction() #=============================================================================== package_get_element_types(core _types) foreach(_type ${_types}) if(EXISTS ${CMAKE_CURRENT_SOURCE_DIR}/${_type}.msh) list(APPEND _meshes ${_type}.msh) register_fem_test(fe_engine_precomputation ${_type}) else() if(NOT ${_type} STREQUAL _point_1) message("The mesh ${_type}.msh is missing, the fe_engine test cannot be activated without it") endif() endif() endforeach() -#register_test(test_interpolate_bernoulli_beam_2 test_interpolate_bernoulli_beam_2.cc) +register_test(test_fe_engine_precomupation_bernoulli_beam_2 PACKAGE core SOURCES test_fe_engine_precomputation_bernoulli_2.cc) #add_mesh(test_fem_circle_1_mesh circle.geo 2 1 OUTPUT circle1.msh) #add_mesh(test_fem_circle_2_mesh circle.geo 2 2 OUTPUT circle2.msh) # Tests for class MeshData macro(register_typed_test test_name type value1 value2) set(target test_${test_name}_${type}) register_test(${target} SOURCES test_${test_name}.cc COMPILE_OPTIONS "TYPE=${type};VALUE1=${value1};VALUE2=${value2}" PACKAGE core ) endmacro() register_typed_test(mesh_data string \"5\" \"10\") register_typed_test(mesh_data UInt 5 10) add_mesh(test_boundary_msh cube.geo 3 1) add_mesh(test_boundary_msh_physical_names cube_physical_names.geo 3 1) register_test(test_mesh_boundary SOURCES test_mesh_boundary.cc DEPENDS test_boundary_msh test_boundary_msh_physical_names PACKAGE core) register_test(test_facet_element_mapping SOURCES test_facet_element_mapping.cc DEPENDS test_boundary_msh_physical_names PACKAGE core) register_gtest_sources( SOURCES test_fe_engine_gauss_integration.cc PACKAGE core ) register_gtest_sources( SOURCES test_gradient.cc PACKAGE core ) register_gtest_sources( SOURCES test_integrate.cc PACKAGE core ) register_gtest_sources( SOURCES test_inverse_map.cc PACKAGE core ) register_gtest_test(test_fe_engine FILES_TO_COPY ${_meshes}) diff --git a/test/test_fe_engine/test_fe_engine_precomputation_bernoulli_2.cc b/test/test_fe_engine/test_fe_engine_precomputation_bernoulli_2.cc new file mode 100644 index 000000000..30644bea7 --- /dev/null +++ b/test/test_fe_engine/test_fe_engine_precomputation_bernoulli_2.cc @@ -0,0 +1,72 @@ +/** + * @file test_fe_engine_precomputation.cc + * + * @author Nicolas Richart + * + * @date creation: Mon Jun 14 2010 + * @date last modification: Mon Jul 13 2015 + * + * @brief test of the fem class + * + * @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 "fe_engine.hh" +#include "shape_structural.hh" +#include "integrator_gauss.hh" +/* -------------------------------------------------------------------------- */ +#include +/* -------------------------------------------------------------------------- */ +using namespace akantu; + +int main(int argc, char *argv[]) { + akantu::initialize(argc, argv); + + // debug::setDebugLevel(dblTest); + + constexpr ElementType type = _bernoulli_beam_2; + UInt dim = ElementClass::getSpatialDimension(); + + Mesh mesh(dim); + Vector node = {0, 0}; + mesh.getNodes().push_back(node); + node = {1, 1}; + mesh.getNodes().push_back(node); + + mesh.addConnectivityType(type); + auto & connectivity = mesh.getConnectivity(type); + Vector elem = {0, 1}; + connectivity.push_back(elem); + + auto fem = + std::make_unique>( + mesh, dim, "test_fem"); + + fem->initShapeFunctions(); + + // std::cout << *fem << std::endl; + + // delete fem; + // finalize(); + + return 0; +}