diff --git a/packages/core.cmake b/packages/core.cmake index 70e2c03fa..0087ce7ab 100644 --- a/packages/core.cmake +++ b/packages/core.cmake @@ -1,492 +1,493 @@ #=============================================================================== # @file core.cmake # # @author Guillaume Anciaux # @author Nicolas Richart # # @date creation: Mon Nov 21 2011 # @date last modification: Mon Jan 18 2016 # # @brief package description for core # # @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 . # #=============================================================================== package_declare(core NOT_OPTIONAL DESCRIPTION "core package for Akantu") package_declare_sources(core common/aka_array.cc common/aka_array.hh common/aka_array_tmpl.hh common/aka_blas_lapack.hh common/aka_circular_array.hh common/aka_circular_array_inline_impl.cc common/aka_common.cc common/aka_common.hh common/aka_common_inline_impl.cc common/aka_csr.hh common/aka_element_classes_info_inline_impl.cc common/aka_error.cc common/aka_error.hh common/aka_event_handler_manager.hh common/aka_extern.cc common/aka_fwd.hh common/aka_grid_dynamic.hh common/aka_math.cc common/aka_math.hh common/aka_math_tmpl.hh common/aka_memory.cc common/aka_memory.hh common/aka_memory_inline_impl.cc common/aka_random_generator.hh common/aka_safe_enum.hh common/aka_static_memory.cc common/aka_static_memory.hh common/aka_static_memory_inline_impl.cc common/aka_static_memory_tmpl.hh common/aka_typelist.hh common/aka_types.hh common/aka_visitor.hh common/aka_voigthelper.hh common/aka_voigthelper.cc fe_engine/element_class.cc fe_engine/element_class.hh fe_engine/element_class_tmpl.hh fe_engine/element_classes/element_class_hexahedron_8_inline_impl.cc fe_engine/element_classes/element_class_hexahedron_20_inline_impl.cc fe_engine/element_classes/element_class_pentahedron_6_inline_impl.cc fe_engine/element_classes/element_class_pentahedron_15_inline_impl.cc fe_engine/element_classes/element_class_point_1_inline_impl.cc fe_engine/element_classes/element_class_quadrangle_4_inline_impl.cc fe_engine/element_classes/element_class_quadrangle_8_inline_impl.cc fe_engine/element_classes/element_class_segment_2_inline_impl.cc fe_engine/element_classes/element_class_segment_3_inline_impl.cc fe_engine/element_classes/element_class_tetrahedron_10_inline_impl.cc fe_engine/element_classes/element_class_tetrahedron_4_inline_impl.cc fe_engine/element_classes/element_class_triangle_3_inline_impl.cc fe_engine/element_classes/element_class_triangle_6_inline_impl.cc fe_engine/fe_engine.cc fe_engine/fe_engine.hh fe_engine/fe_engine_inline_impl.cc fe_engine/fe_engine_template.hh fe_engine/fe_engine_template_tmpl.hh fe_engine/geometrical_element.cc - fe_engine/integration_element.cc + fe_engine/gauss_integration.cc + fe_engine/gauss_integration_tmpl.hh fe_engine/integrator.hh fe_engine/integrator_gauss.hh fe_engine/integrator_gauss_inline_impl.cc fe_engine/interpolation_element.cc fe_engine/interpolation_element_tmpl.hh fe_engine/integration_point.hh fe_engine/shape_functions.hh fe_engine/shape_functions_inline_impl.cc fe_engine/shape_lagrange.cc fe_engine/shape_lagrange.hh fe_engine/shape_lagrange_inline_impl.cc fe_engine/shape_linked.cc fe_engine/shape_linked.hh fe_engine/shape_linked_inline_impl.cc fe_engine/element.hh io/dumper/dumpable.hh io/dumper/dumpable.cc io/dumper/dumpable_dummy.hh io/dumper/dumpable_inline_impl.hh io/dumper/dumper_field.hh io/dumper/dumper_material_padders.hh io/dumper/dumper_filtered_connectivity.hh io/dumper/dumper_element_partition.hh io/mesh_io.cc io/mesh_io.hh io/mesh_io/mesh_io_abaqus.cc io/mesh_io/mesh_io_abaqus.hh io/mesh_io/mesh_io_diana.cc io/mesh_io/mesh_io_diana.hh io/mesh_io/mesh_io_msh.cc io/mesh_io/mesh_io_msh.hh io/model_io.cc io/model_io.hh io/parser/algebraic_parser.hh io/parser/input_file_parser.hh io/parser/parsable.cc io/parser/parsable.hh io/parser/parsable_tmpl.hh io/parser/parser.cc io/parser/parser_real.cc io/parser/parser_random.cc io/parser/parser_types.cc io/parser/parser_input_files.cc io/parser/parser.hh io/parser/parser_tmpl.hh io/parser/parser_grammar_tmpl.hh io/parser/cppargparse/cppargparse.hh io/parser/cppargparse/cppargparse.cc io/parser/cppargparse/cppargparse_tmpl.hh mesh/element_group.cc mesh/element_group.hh mesh/element_group_inline_impl.cc mesh/element_type_map.hh mesh/element_type_map_tmpl.hh mesh/element_type_map_filter.hh mesh/group_manager.cc mesh/group_manager.hh mesh/group_manager_inline_impl.cc mesh/mesh.cc mesh/mesh.hh mesh/mesh_accessor.hh mesh/mesh_events.hh mesh/mesh_filter.hh mesh/mesh_data.cc mesh/mesh_data.hh mesh/mesh_data_tmpl.hh mesh/mesh_inline_impl.cc mesh/node_group.cc mesh/node_group.hh mesh/node_group_inline_impl.cc mesh_utils/mesh_partition.cc mesh_utils/mesh_partition.hh mesh_utils/mesh_partition/mesh_partition_mesh_data.cc mesh_utils/mesh_partition/mesh_partition_mesh_data.hh mesh_utils/mesh_partition/mesh_partition_scotch.hh mesh_utils/mesh_utils_pbc.cc mesh_utils/mesh_utils.cc mesh_utils/mesh_utils.hh mesh_utils/mesh_utils_inline_impl.cc mesh_utils/global_ids_updater.hh mesh_utils/global_ids_updater.cc mesh_utils/global_ids_updater_inline_impl.cc model/boundary_condition.hh model/boundary_condition_functor.hh model/boundary_condition_functor_inline_impl.cc model/boundary_condition_tmpl.hh model/integration_scheme/generalized_trapezoidal.hh model/integration_scheme/generalized_trapezoidal_inline_impl.cc model/integration_scheme/integration_scheme_1st_order.hh model/integration_scheme/integration_scheme_2nd_order.hh model/integration_scheme/newmark-beta.hh model/integration_scheme/newmark-beta_inline_impl.cc model/model.cc model/model.hh model/model_inline_impl.cc model/solid_mechanics/material.cc model/solid_mechanics/material.hh model/solid_mechanics/material_inline_impl.cc model/solid_mechanics/material_selector.hh model/solid_mechanics/material_selector_tmpl.hh model/solid_mechanics/materials/internal_field.hh model/solid_mechanics/materials/internal_field_tmpl.hh model/solid_mechanics/materials/random_internal_field.hh model/solid_mechanics/materials/random_internal_field_tmpl.hh model/solid_mechanics/solid_mechanics_model.cc model/solid_mechanics/solid_mechanics_model.hh model/solid_mechanics/solid_mechanics_model_inline_impl.cc model/solid_mechanics/solid_mechanics_model_mass.cc model/solid_mechanics/solid_mechanics_model_material.cc model/solid_mechanics/solid_mechanics_model_tmpl.hh model/solid_mechanics/solid_mechanics_model_event_handler.hh model/solid_mechanics/materials/plane_stress_toolbox.hh model/solid_mechanics/materials/plane_stress_toolbox_tmpl.hh model/solid_mechanics/materials/material_core_includes.hh model/solid_mechanics/materials/material_elastic.cc model/solid_mechanics/materials/material_elastic.hh model/solid_mechanics/materials/material_elastic_inline_impl.cc model/solid_mechanics/materials/material_thermal.cc model/solid_mechanics/materials/material_thermal.hh model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc model/solid_mechanics/materials/material_elastic_linear_anisotropic.hh model/solid_mechanics/materials/material_elastic_orthotropic.cc model/solid_mechanics/materials/material_elastic_orthotropic.hh model/solid_mechanics/materials/material_damage/material_damage.hh model/solid_mechanics/materials/material_damage/material_damage_tmpl.hh model/solid_mechanics/materials/material_damage/material_marigo.cc model/solid_mechanics/materials/material_damage/material_marigo.hh model/solid_mechanics/materials/material_damage/material_marigo_inline_impl.cc model/solid_mechanics/materials/material_damage/material_mazars.cc model/solid_mechanics/materials/material_damage/material_mazars.hh model/solid_mechanics/materials/material_damage/material_mazars_inline_impl.cc model/solid_mechanics/materials/material_finite_deformation/material_neohookean.cc model/solid_mechanics/materials/material_finite_deformation/material_neohookean.hh model/solid_mechanics/materials/material_finite_deformation/material_neohookean_inline_impl.cc model/solid_mechanics/materials/material_plastic/material_plastic.cc model/solid_mechanics/materials/material_plastic/material_plastic.hh model/solid_mechanics/materials/material_plastic/material_plastic_inline_impl.cc model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening.cc model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening.hh model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening_inline_impl.cc model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.cc model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.hh model/common/neighborhood_base.hh model/common/neighborhood_base.cc model/common/neighborhood_base_inline_impl.cc model/common/neighborhoods_criterion_evaluation/neighborhood_max_criterion.hh model/common/neighborhoods_criterion_evaluation/neighborhood_max_criterion.cc model/common/neighborhoods_criterion_evaluation/neighborhood_max_criterion_inline_impl.cc solver/solver.cc solver/solver.hh solver/solver_inline_impl.cc solver/sparse_matrix.cc solver/sparse_matrix.hh solver/sparse_matrix_inline_impl.cc solver/static_solver.hh solver/static_solver.cc synchronizer/communication_buffer.hh synchronizer/communication_buffer_inline_impl.cc synchronizer/data_accessor.cc synchronizer/data_accessor.hh synchronizer/data_accessor_inline_impl.cc synchronizer/distributed_synchronizer.cc synchronizer/distributed_synchronizer.hh synchronizer/distributed_synchronizer_tmpl.hh synchronizer/dof_synchronizer.cc synchronizer/dof_synchronizer.hh synchronizer/dof_synchronizer_inline_impl.cc synchronizer/filtered_synchronizer.cc synchronizer/filtered_synchronizer.hh synchronizer/pbc_synchronizer.cc synchronizer/pbc_synchronizer.hh synchronizer/real_static_communicator.hh synchronizer/static_communicator.cc synchronizer/static_communicator.hh synchronizer/static_communicator_dummy.hh synchronizer/static_communicator_inline_impl.hh synchronizer/synchronizer.cc synchronizer/synchronizer.hh synchronizer/synchronizer_registry.cc synchronizer/synchronizer_registry.hh synchronizer/grid_synchronizer.cc synchronizer/grid_synchronizer.hh ) package_declare_elements(core ELEMENT_TYPES _point_1 _segment_2 _segment_3 _triangle_3 _triangle_6 _quadrangle_4 _quadrangle_8 _tetrahedron_4 _tetrahedron_10 _pentahedron_6 _pentahedron_15 _hexahedron_8 _hexahedron_20 KIND regular GEOMETRICAL_TYPES _gt_point _gt_segment_2 _gt_segment_3 _gt_triangle_3 _gt_triangle_6 _gt_quadrangle_4 _gt_quadrangle_8 _gt_tetrahedron_4 _gt_tetrahedron_10 _gt_hexahedron_8 _gt_hexahedron_20 _gt_pentahedron_6 _gt_pentahedron_15 INTERPOLATION_TYPES _itp_lagrange_point_1 _itp_lagrange_segment_2 _itp_lagrange_segment_3 _itp_lagrange_triangle_3 _itp_lagrange_triangle_6 _itp_lagrange_quadrangle_4 _itp_serendip_quadrangle_8 _itp_lagrange_tetrahedron_4 _itp_lagrange_tetrahedron_10 _itp_lagrange_hexahedron_8 _itp_serendip_hexahedron_20 _itp_lagrange_pentahedron_6 _itp_lagrange_pentahedron_15 GEOMETRICAL_SHAPES _gst_point _gst_triangle _gst_square _gst_prism GAUSS_INTEGRATION_TYPES _git_point _git_segment _git_triangle _git_tetrahedron _git_pentahedron INTERPOLATION_KIND _itk_lagrangian FE_ENGINE_LISTS gradient_on_integration_points interpolate_on_integration_points interpolate compute_normals_on_integration_points inverse_map contains compute_shapes compute_shapes_derivatives get_shapes_derivatives ) package_declare_material_infos(core LIST AKANTU_CORE_MATERIAL_LIST INCLUDE material_core_includes.hh ) package_declare_documentation_files(core manual.sty manual.cls manual.tex manual-macros.sty manual-titlepages.tex manual-introduction.tex manual-gettingstarted.tex manual-io.tex manual-feengine.tex manual-solidmechanicsmodel.tex manual-constitutive-laws.tex manual-lumping.tex manual-elements.tex manual-appendix-elements.tex manual-appendix-materials.tex manual-appendix-packages.tex manual-backmatter.tex manual-bibliography.bib manual-bibliographystyle.bst figures/bc_and_ic_example.pdf figures/boundary.pdf figures/boundary.svg figures/dirichlet.pdf figures/dirichlet.svg figures/doc_wheel.pdf figures/doc_wheel.svg figures/dynamic_analysis.png figures/explicit_dynamic.pdf figures/explicit_dynamic.svg figures/static.pdf figures/static.svg figures/hooke_law.pdf figures/hot-point-1.png figures/hot-point-2.png figures/implicit_dynamic.pdf figures/implicit_dynamic.svg figures/insertion.pdf figures/interpolate.pdf figures/interpolate.svg figures/problemDomain.pdf_tex figures/problemDomain.pdf figures/static_analysis.png figures/stress_strain_el.pdf figures/tangent.pdf figures/tangent.svg figures/vectors.pdf figures/vectors.svg figures/stress_strain_neo.pdf figures/visco_elastic_law.pdf figures/isotropic_hardening_plasticity.pdf figures/stress_strain_visco.pdf figures/elements/hexahedron_8.pdf figures/elements/hexahedron_8.svg figures/elements/quadrangle_4.pdf figures/elements/quadrangle_4.svg figures/elements/quadrangle_8.pdf figures/elements/quadrangle_8.svg figures/elements/segment_2.pdf figures/elements/segment_2.svg figures/elements/segment_3.pdf figures/elements/segment_3.svg figures/elements/tetrahedron_10.pdf figures/elements/tetrahedron_10.svg figures/elements/tetrahedron_4.pdf figures/elements/tetrahedron_4.svg figures/elements/triangle_3.pdf figures/elements/triangle_3.svg figures/elements/triangle_6.pdf figures/elements/triangle_6.svg figures/elements/xtemp.pdf ) package_declare_documentation(core "This package is the core engine of \\akantu. It depends on:" "\\begin{itemize}" "\\item A C++ compiler (\\href{http://gcc.gnu.org/}{GCC} >= 4, or \\href{https://software.intel.com/en-us/intel-compilers}{Intel})." "\\item The cross-platform, open-source \\href{http://www.cmake.org/}{CMake} build system." "\\item The \\href{http://www.boost.org/}{Boost} C++ portable libraries." "\\item The \\href{http://www.zlib.net/}{zlib} compression library." "\\end{itemize}" "" "Under Ubuntu (14.04 LTS) the installation can be performed using the commands:" "\\begin{command}" " > sudo apt-get install cmake libboost-dev zlib1g-dev g++" "\\end{command}" "" "Under Mac OS X the installation requires the following steps:" "\\begin{itemize}" "\\item Install Xcode" "\\item Install the command line tools." "\\item Install the MacPorts project which allows to automatically" "download and install opensource packages." "\\end{itemize}" "Then the following commands should be typed in a terminal:" "\\begin{command}" " > sudo port install cmake gcc48 boost" "\\end{command}" ) find_program(READLINK_COMMAND readlink) find_program(ADDR2LINE_COMMAND addr2line) find_program(PATCH_COMMAND patch) mark_as_advanced(READLINK_COMMAND) mark_as_advanced(ADDR2LINE_COMMAND) include(CheckFunctionExists) check_function_exists(clock_gettime _clock_gettime) include(CheckCXXSymbolExists) check_cxx_symbol_exists(strdup cstring AKANTU_HAS_STRDUP) if(NOT _clock_gettime) set(AKANTU_USE_OBSOLETE_GETTIMEOFDAY ON CACHE INTERNAL "" FORCE) else() set(AKANTU_USE_OBSOLETE_GETTIMEOFDAY OFF CACHE INTERNAL "" FORCE) endif() package_declare_extra_files_to_package(core SOURCES common/aka_element_classes_info.hh.in common/aka_config.hh.in model/solid_mechanics/material_list.hh.in ) diff --git a/src/fe_engine/element_class_tmpl.hh b/src/fe_engine/element_class_tmpl.hh index 0203d9291..272116bca 100644 --- a/src/fe_engine/element_class_tmpl.hh +++ b/src/fe_engine/element_class_tmpl.hh @@ -1,676 +1,512 @@ /** * @file element_class_tmpl.hh * * @author Aurelia Isabel Cuba Ramos * @author Thomas Menouillard * @author Nicolas Richart * * @date creation: Thu Feb 21 2013 * @date last modification: Thu Oct 22 2015 * * @brief Implementation of the inline templated function of the element class * descriptions * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifdef AKANTU_CORE_CXX11 #include #endif /* -------------------------------------------------------------------------- */ -__BEGIN_AKANTU__ - -/* -------------------------------------------------------------------------- */ -/* GaussIntegrationElement */ -/* -------------------------------------------------------------------------- */ -namespace _aka_gauss_helpers { - template struct GaussIntegrationNbPoints { - static const UInt nb_points = -1; - }; - - template struct GaussIntegrationNbPoints<_git_point, n> { - static const UInt nb_points = 1; - }; - - template struct GaussIntegrationNbPoints<_git_segment, n> { - static const UInt nb_points = (n + 1) / 2 + ((n + 1) % 2 ? 1 : 0); - }; - -#define DECLARE_GAUSS_NB_POINTS(type, order, points) \ - template <> struct GaussIntegrationNbPoints { \ - static const UInt nb_points = points; \ - } - - DECLARE_GAUSS_NB_POINTS(_git_triangle, 1, 1); - DECLARE_GAUSS_NB_POINTS(_git_triangle, 2, 3); - DECLARE_GAUSS_NB_POINTS(_git_triangle, 3, 4); - DECLARE_GAUSS_NB_POINTS(_git_triangle, 4, 6); - DECLARE_GAUSS_NB_POINTS(_git_tetrahedron, 1, 1); - DECLARE_GAUSS_NB_POINTS(_git_tetrahedron, 2, 4); - DECLARE_GAUSS_NB_POINTS(_git_tetrahedron, 3, 5); - DECLARE_GAUSS_NB_POINTS(_git_tetrahedron, 5, 15); - DECLARE_GAUSS_NB_POINTS(_git_pentahedron, 2, 6); // order 3 in x, order 2 in y and z - DECLARE_GAUSS_NB_POINTS(_git_pentahedron, 3, 8); // order 3 in x, order 3 in y and z - DECLARE_GAUSS_NB_POINTS(_git_pentahedron, 5, 21); // order 5 in x, order 5 in y and z - - template - struct GaussIntegrationNbPointsHelper { - static const UInt pnp = GaussIntegrationNbPoints::nb_points; - static const UInt nb_points = - ((pnp == UInt(-1)) ? - GaussIntegrationNbPointsHelper 2 * on)>::nb_points - : pnp); - }; - - template - struct GaussIntegrationNbPointsHelper { - static const UInt nb_points = 0; - }; - - template - struct GaussIntegrationTypeMissingData { - typedef void value; - }; - - template - struct GaussIntegrationTypeMissingData { - // If you fall here it means you are missing GaussIntegrationTypeData for - // the GaussIntegrationType 'type' with polynomials of order 'n' - }; - - /* ------------------------------------------------------------------------ */ - template - struct GaussIntegrationTypeDataHelper { - typedef GaussIntegrationNbPointsHelper git_np; - typedef GaussIntegrationTypeData git_data; - typedef typename GaussIntegrationTypeMissingData< - type, n, git_np::nb_points>::value git_missing; - - static UInt getNbQuadraturePoints() { return git_np::nb_points; } - - static const Matrix getQuadraturePoints() { - return Matrix(git_data::quad_positions, dimension, git_np::nb_points); - } - - static const Vector getWeights() { - return Vector(git_data::quad_weights, git_np::nb_points); - } - }; - - /* -------------------------------------------------------------------------- */ - template - struct GaussIntegrationTypeDataHelper<_git_segment, dimension, dp> { - typedef GaussIntegrationNbPointsHelper<_git_segment, dp> git_np; - typedef GaussIntegrationTypeData<_git_segment, git_np::nb_points> git_data; - typedef typename GaussIntegrationTypeMissingData< - _git_segment, dp, git_np::nb_points>::value git_missing; - - static UInt getNbQuadraturePoints() { - return Math::pow(git_np::nb_points); - } - - static const Matrix getQuadraturePoints() { - UInt tot_nquad = getNbQuadraturePoints(); - UInt nquad = git_np::nb_points; - - Matrix quads(dimension, tot_nquad); - Vector pos(git_data::quad_positions, nquad); - - UInt offset = 1; - for (UInt d = 0; d < dimension; ++d) { - for (UInt n = 0, q = 0; n < tot_nquad; ++n, q += offset) { - UInt rq = q % tot_nquad + q / tot_nquad; - quads(d, rq) = pos(n % nquad); - } - offset *= nquad; - } - return quads; - } +#include "gauss_integration_tmpl.hh" - static const Vector getWeights() { - UInt tot_nquad = getNbQuadraturePoints(); - UInt nquad = git_np::nb_points; - - Vector quads_weights(tot_nquad, 1.); - Vector weights(git_data::quad_weights, nquad); - - UInt offset = 1; - for (UInt d = 0; d < dimension; ++d) { - for (UInt n = 0, q = 0; n < tot_nquad; ++n, q += offset) { - UInt rq = q % tot_nquad + q / tot_nquad; - quads_weights(rq) *= weights(n % nquad); - } - offset *= nquad; - } - return quads_weights; - } - }; -} - -template -const Matrix -GaussIntegrationElement::getQuadraturePoints() { - const InterpolationType itp_type = - ElementClassProperty::interpolation_type; - typedef InterpolationPorperty interpolation_property; - typedef _aka_gauss_helpers::GaussIntegrationTypeDataHelper< - ElementClassProperty::gauss_integration_type, - interpolation_property::natural_space_dimension, n> data_helper; - Matrix tmp(data_helper::getQuadraturePoints()); - return tmp; -} - -/* -------------------------------------------------------------------------- */ -template -const Vector GaussIntegrationElement::getWeights() { - const InterpolationType itp_type = - ElementClassProperty::interpolation_type; - typedef InterpolationPorperty interpolation_property; - typedef _aka_gauss_helpers::GaussIntegrationTypeDataHelper< - ElementClassProperty::gauss_integration_type, - interpolation_property::natural_space_dimension, n> data_helper; - Vector tmp(data_helper::getWeights()); - return tmp; -} - -/* -------------------------------------------------------------------------- */ -template -UInt GaussIntegrationElement::getNbQuadraturePoints() { - const InterpolationType itp_type = - ElementClassProperty::interpolation_type; - typedef InterpolationPorperty interpolation_property; - typedef _aka_gauss_helpers::GaussIntegrationTypeDataHelper< - ElementClassProperty::gauss_integration_type, - interpolation_property::natural_space_dimension, n> data_helper; - return data_helper::getNbQuadraturePoints(); -} +__BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* GeometricalElement */ /* -------------------------------------------------------------------------- */ template inline const MatrixProxy GeometricalElement::getFacetLocalConnectivityPerElement(UInt t) { return MatrixProxy(facet_connectivity[t], nb_facets[t], nb_nodes_per_facet[t]); } /* -------------------------------------------------------------------------- */ template inline UInt GeometricalElement::getNbFacetsPerElement() { UInt total_nb_facets = 0; for (UInt n = 0; n < nb_facet_types; ++n) { total_nb_facets += nb_facets[n]; } return total_nb_facets; } /* -------------------------------------------------------------------------- */ template inline UInt GeometricalElement::getNbFacetsPerElement(UInt t) { return nb_facets[t]; } /* -------------------------------------------------------------------------- */ template template inline bool GeometricalElement::contains( const vector_type & coords) { return GeometricalShapeContains::contains(coords); } /* -------------------------------------------------------------------------- */ template <> template inline bool GeometricalShapeContains<_gst_point>::contains(const vector_type & coords) { return (coords(0) < std::numeric_limits::epsilon()); } /* -------------------------------------------------------------------------- */ template <> template inline bool GeometricalShapeContains<_gst_square>::contains(const vector_type & coords) { bool in = true; for (UInt i = 0; i < coords.size() && in; ++i) in &= ((coords(i) >= -(1. + std::numeric_limits::epsilon())) && (coords(i) <= (1. + std::numeric_limits::epsilon()))); return in; } /* -------------------------------------------------------------------------- */ template <> template inline bool GeometricalShapeContains<_gst_triangle>::contains(const vector_type & coords) { bool in = true; Real sum = 0; for (UInt i = 0; (i < coords.size()) && in; ++i) { in &= ((coords(i) >= -(Math::getTolerance())) && (coords(i) <= (1. + Math::getTolerance()))); sum += coords(i); } if (in) return (in && (sum <= (1. + Math::getTolerance()))); return in; } /* -------------------------------------------------------------------------- */ template <> template inline bool GeometricalShapeContains<_gst_prism>::contains(const vector_type & coords) { bool in = ((coords(0) >= -1.) && (coords(0) <= 1.)); // x in segment [-1, 1] // y and z in triangle in &= ((coords(1) >= 0) && (coords(1) <= 1.)); in &= ((coords(2) >= 0) && (coords(2) <= 1.)); Real sum = coords(1) + coords(2); return (in && (sum <= 1)); } /* -------------------------------------------------------------------------- */ /* InterpolationElement */ /* -------------------------------------------------------------------------- */ template inline void InterpolationElement::computeShapes( const Matrix & natural_coord, Matrix & N) { UInt nb_points = natural_coord.cols(); for (UInt p = 0; p < nb_points; ++p) { Vector Np(N(p)); Vector ncoord_p(natural_coord(p)); computeShapes(ncoord_p, Np); } } /* -------------------------------------------------------------------------- */ template inline void InterpolationElement::computeDNDS( const Matrix & natural_coord, Tensor3 & dnds) { UInt nb_points = natural_coord.cols(); for (UInt p = 0; p < nb_points; ++p) { Matrix dnds_p(dnds(p)); Vector ncoord_p(natural_coord(p)); computeDNDS(ncoord_p, dnds_p); } } /* -------------------------------------------------------------------------- */ /** * interpolate on a point a field for which values are given on the * node of the element using the shape functions at this interpolation point * * @param nodal_values values of the function per node @f$ f_{ij} = f_{n_i j} *@f$ so it should be a matrix of size nb_nodes_per_element @f$\times@f$ *nb_degree_of_freedom * @param shapes value of shape functions at the interpolation point * @param interpolated interpolated value of f @f$ f_j(\xi) = \sum_i f_{n_i j} *N_i @f$ */ template inline void InterpolationElement::interpolate( const Matrix & nodal_values, const Vector & shapes, Vector & interpolated) { Matrix interpm(interpolated.storage(), nodal_values.rows(), 1); Matrix shapesm( shapes.storage(), InterpolationPorperty::nb_nodes_per_element, 1); interpm.mul(nodal_values, shapesm); } /* -------------------------------------------------------------------------- */ /** * interpolate on several points a field for which values are given on the * node of the element using the shape functions at the interpolation point * * @param nodal_values values of the function per node @f$ f_{ij} = f_{n_i j} *@f$ so it should be a matrix of size nb_nodes_per_element @f$\times@f$ *nb_degree_of_freedom * @param shapes value of shape functions at the interpolation point * @param interpolated interpolated values of f @f$ f_j(\xi) = \sum_i f_{n_i j} *N_i @f$ */ template inline void InterpolationElement::interpolate( const Matrix & nodal_values, const Matrix & shapes, Matrix & interpolated) { UInt nb_points = shapes.cols(); for (UInt p = 0; p < nb_points; ++p) { Vector Np(shapes(p)); Vector interpolated_p(interpolated(p)); interpolate(nodal_values, Np, interpolated_p); } } /* -------------------------------------------------------------------------- */ /** * interpolate the field on a point given in natural coordinates the field which * values are given on the node of the element * * @param natural_coords natural coordinates of point where to interpolate \xi * @param nodal_values values of the function per node @f$ f_{ij} = f_{n_i j} *@f$ so it should be a matrix of size nb_nodes_per_element @f$\times@f$ *nb_degree_of_freedom * @param interpolated interpolated value of f @f$ f_j(\xi) = \sum_i f_{n_i j} *N_i @f$ */ template inline void InterpolationElement::interpolateOnNaturalCoordinates( const Vector & natural_coords, const Matrix & nodal_values, Vector & interpolated) { Vector shapes( InterpolationPorperty::nb_nodes_per_element); computeShapes(natural_coords, shapes); interpolate(nodal_values, shapes, interpolated); } /* -------------------------------------------------------------------------- */ /// @f$ gradient_{ij} = \frac{\partial f_j}{\partial s_i} = \sum_k /// \frac{\partial N_k}{\partial s_i}f_{j n_k} @f$ template inline void InterpolationElement::gradientOnNaturalCoordinates( const Vector & natural_coords, const Matrix & f, Matrix & gradient) { Matrix dnds( InterpolationPorperty::natural_space_dimension, InterpolationPorperty::nb_nodes_per_element); computeDNDS(natural_coords, dnds); gradient.mul(f, dnds); } /* -------------------------------------------------------------------------- */ /* ElementClass */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeJMat(const Tensor3 & dnds, const Matrix & node_coords, Tensor3 & J) { UInt nb_points = dnds.size(2); for (UInt p = 0; p < nb_points; ++p) { Matrix J_p(J(p)); Matrix dnds_p(dnds(p)); computeJMat(dnds_p, node_coords, J_p); } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeJMat(const Matrix & dnds, const Matrix & node_coords, Matrix & J) { /// @f$ J = dxds = dnds * x @f$ J.mul(dnds, node_coords); } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeJacobian(const Matrix & natural_coords, const Matrix & node_coords, Vector & jacobians) { UInt nb_points = natural_coords.cols(); Matrix dnds(interpolation_property::natural_space_dimension, interpolation_property::nb_nodes_per_element); Matrix J(natural_coords.rows(), node_coords.rows()); for (UInt p = 0; p < nb_points; ++p) { Vector ncoord_p(natural_coords(p)); interpolation_element::computeDNDS(ncoord_p, dnds); computeJMat(dnds, node_coords, J); computeJacobian(J, jacobians(p)); } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeJacobian(const Tensor3 & J, Vector & jacobians) { UInt nb_points = J.size(2); for (UInt p = 0; p < nb_points; ++p) { computeJacobian(J(p), jacobians(p)); } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeJacobian(const Matrix & J, Real & jacobians) { if (J.rows() == J.cols()) { jacobians = Math::det(J.storage()); } else { interpolation_element::computeSpecialJacobian(J, jacobians); } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapeDerivatives(const Tensor3 & J, const Tensor3 & dnds, Tensor3 & shape_deriv) { UInt nb_points = J.size(2); for (UInt p = 0; p < nb_points; ++p) { Matrix shape_deriv_p(shape_deriv(p)); computeShapeDerivatives(J(p), dnds(p), shape_deriv_p); } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapeDerivatives(const Matrix & J, const Matrix & dnds, Matrix & shape_deriv) { Matrix inv_J(J.rows(), J.cols()); Math::inv(J.storage(), inv_J.storage()); shape_deriv.mul(inv_J, dnds); } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeNormalsOnNaturalCoordinates( const Matrix & coord, Matrix & f, Matrix & normals) { UInt dimension = normals.rows(); UInt nb_points = coord.cols(); AKANTU_DEBUG_ASSERT((dimension - 1) == interpolation_property::natural_space_dimension, "cannot extract a normal because of dimension mismatch " << dimension - 1 << " " << interpolation_property::natural_space_dimension); Matrix J(dimension, interpolation_property::natural_space_dimension); for (UInt p = 0; p < nb_points; ++p) { interpolation_element::gradientOnNaturalCoordinates(coord(p), f, J); if (dimension == 2) { Math::normal2(J.storage(), normals(p).storage()); } if (dimension == 3) { Math::normal3(J(0).storage(), J(1).storage(), normals(p).storage()); } } } /* ------------------------------------------------------------------------- */ /** * In the non linear cases we need to iterate to find the natural coordinates *@f$\xi@f$ * provided real coordinates @f$x@f$. * * We want to solve: @f$ x- \phi(\xi) = 0@f$ with @f$\phi(\xi) = \sum_I N_I(\xi) *x_I@f$ * the mapping function which uses the nodal coordinates @f$x_I@f$. * * To that end we use the Newton method and the following series: * * @f$ \frac{\partial \phi(x_k)}{\partial \xi} \left( \xi_{k+1} - \xi_k \right) *= x - \phi(x_k)@f$ * * When we consider elements embedded in a dimension higher than them (2D *triangle in a 3D space for example) * @f$ J = \frac{\partial \phi(\xi_k)}{\partial \xi}@f$ is of dimension *@f$dim_{space} \times dim_{elem}@f$ which * is not invertible in most cases. Rather we can solve the problem: * * @f$ J^T J \left( \xi_{k+1} - \xi_k \right) = J^T \left( x - \phi(\xi_k) *\right) @f$ * * So that * * @f$ d\xi = \xi_{k+1} - \xi_k = (J^T J)^{-1} J^T \left( x - \phi(\xi_k) *\right) @f$ * * So that if the series converges we have: * * @f$ 0 = J^T \left( \phi(\xi_\infty) - x \right) @f$ * * And we see that this is ill-posed only if @f$ J^T x = 0@f$ which means that *the vector provided * is normal to any tangent which means it is outside of the element itself. * * @param real_coords: the real coordinates the natural coordinates are sought *for * @param node_coords: the coordinates of the nodes forming the element * @param natural_coords: output->the sought natural coordinates * @param spatial_dimension: spatial dimension of the problem * **/ template inline void ElementClass::inverseMap( const Vector & real_coords, const Matrix & node_coords, Vector & natural_coords, Real tolerance) { UInt spatial_dimension = real_coords.size(); UInt dimension = natural_coords.size(); // matrix copy of the real_coords Matrix mreal_coords(real_coords.storage(), spatial_dimension, 1); // initial guess // Matrix natural_guess(natural_coords.storage(), dimension, 1); natural_coords.clear(); // real space coordinates provided by initial guess Matrix physical_guess(dimension, 1); // objective function f = real_coords - physical_guess Matrix f(dimension, 1); // dnds computed on the natural_guess // Matrix dnds(interpolation_element::nb_nodes_per_element, // spatial_dimension); // J Jacobian matrix computed on the natural_guess Matrix J(spatial_dimension, dimension); // G = J^t * J Matrix G(spatial_dimension, spatial_dimension); // Ginv = G^{-1} Matrix Ginv(spatial_dimension, spatial_dimension); // J = Ginv * J^t Matrix F(spatial_dimension, dimension); // dxi = \xi_{k+1} - \xi in the iterative process Matrix dxi(spatial_dimension, 1); /* --------------------------- */ /* init before iteration loop */ /* --------------------------- */ // do interpolation Vector physical_guess_v(physical_guess.storage(), dimension); interpolation_element::interpolateOnNaturalCoordinates( natural_coords, node_coords, physical_guess_v); // compute initial objective function value f = real_coords - physical_guess f = mreal_coords; f -= physical_guess; // compute initial error Real inverse_map_error = f.norm(); /* --------------------------- */ /* iteration loop */ /* --------------------------- */ while (tolerance < inverse_map_error) { // compute J^t interpolation_element::gradientOnNaturalCoordinates(natural_coords, node_coords, J); // compute G G.mul(J, J); // inverse G Ginv.inverse(G); // compute F F.mul(Ginv, J); // compute increment dxi.mul(F, f); // update our guess natural_coords += Vector(dxi(0)); // interpolate interpolation_element::interpolateOnNaturalCoordinates( natural_coords, node_coords, physical_guess_v); // compute error f = mreal_coords; f -= physical_guess; inverse_map_error = f.norm(); } // memcpy(natural_coords.storage(), natural_guess.storage(), sizeof(Real) * // natural_coords.size()); } /* -------------------------------------------------------------------------- */ template inline void ElementClass::inverseMap( const Matrix & real_coords, const Matrix & node_coords, Matrix & natural_coords, Real tolerance) { UInt nb_points = real_coords.cols(); for (UInt p = 0; p < nb_points; ++p) { Vector X(real_coords(p)); Vector ncoord_p(natural_coords(p)); inverseMap(X, node_coords, ncoord_p, tolerance); } } __END_AKANTU__ diff --git a/src/fe_engine/integration_element.cc b/src/fe_engine/gauss_integration.cc similarity index 65% rename from src/fe_engine/integration_element.cc rename to src/fe_engine/gauss_integration.cc index 28e5eb703..5c3571d26 100644 --- a/src/fe_engine/integration_element.cc +++ b/src/fe_engine/gauss_integration.cc @@ -1,218 +1,173 @@ /** * @file integration_element.cc * * @author Mauro Corrado * @author Sébastien Hartmann * @author Thomas Menouillard * @author Nicolas Richart * * @date creation: Thu Feb 21 2013 * @date last modification: Thu Nov 05 2015 * * @brief Definition of the integration constants, some of the value are taken * from r3.01.01 doc from Code Aster * * @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 "aka_common.hh" #include "element_class.hh" using std::sqrt; __BEGIN_AKANTU__ +/* clang-format off */ /* -------------------------------------------------------------------------- */ /* Points */ /* -------------------------------------------------------------------------- */ template<> Real GaussIntegrationTypeData<_git_point, 1>::quad_positions[] = {0}; template<> Real GaussIntegrationTypeData<_git_point, 1>::quad_weights[] = {1.}; /* -------------------------------------------------------------------------- */ /* Segments */ /* -------------------------------------------------------------------------- */ template<> Real GaussIntegrationTypeData<_git_segment, 1>::quad_positions[] = {0.}; template<> Real GaussIntegrationTypeData<_git_segment, 1>::quad_weights[] = {2.}; /* -------------------------------------------------------------------------- */ template<> Real GaussIntegrationTypeData<_git_segment, 2>::quad_positions[] = {-1./sqrt(3.), 1./sqrt(3.)}; template<> Real GaussIntegrationTypeData<_git_segment, 2>::quad_weights[] = {1., 1.}; /* -------------------------------------------------------------------------- */ template<> Real GaussIntegrationTypeData<_git_segment, 3>::quad_positions[] = {-sqrt(3./5.), 0., sqrt(3./5.)}; template<> Real GaussIntegrationTypeData<_git_segment, 3>::quad_weights[] = {5./9., 8./9., 5./9.}; /* -------------------------------------------------------------------------- */ template<> Real GaussIntegrationTypeData<_git_segment, 4>::quad_positions[] = {-sqrt((3. + 2.*sqrt(6./5.))/7.), -sqrt((3. - 2.*sqrt(6./5.))/7.), sqrt((3. - 2.*sqrt(6./5.))/7.), sqrt((3. + 2.*sqrt(6./5.))/7.)}; template<> Real GaussIntegrationTypeData<_git_segment, 4>::quad_weights[] = {(18. - sqrt(30.))/36., (18. + sqrt(30.))/36., (18. + sqrt(30.))/36., (18. - sqrt(30.))/36.}; /* -------------------------------------------------------------------------- */ /* Triangles */ /* -------------------------------------------------------------------------- */ template<> Real GaussIntegrationTypeData<_git_triangle, 1>::quad_positions[] = {1./3., 1./3.}; template<> Real GaussIntegrationTypeData<_git_triangle, 1>::quad_weights[] = {1./2.}; /* -------------------------------------------------------------------------- */ template<> Real GaussIntegrationTypeData<_git_triangle, 3>::quad_positions[] = {1./6., 1./6., 2./3., 1./6., 1./6., 2./3.}; template<> Real GaussIntegrationTypeData<_git_triangle, 3>::quad_weights[] = {1./6., 1./6., 1./6.}; /* -------------------------------------------------------------------------- */ template<> Real GaussIntegrationTypeData<_git_triangle, 4>::quad_positions[] = {1./5., 1./5., 3./5., 1./5., 1./5., 3./5., 1./3., 1./3.}; template<> Real GaussIntegrationTypeData<_git_triangle, 4>::quad_weights[] = {25./(24.*4.), 25./(24.*4.), 25./(24.*4.), -27/(24.*4.)}; /* -------------------------------------------------------------------------- */ /// I give up I cannot find a way to get rational values for this so I sadly /// have to type the pre-computed values at a given precision static const Real tri_6_a = 0.445948490915965; static const Real tri_6_b = 0.091576213509771; static const Real tri_6_w1 = 0.054975871827661; static const Real tri_6_w2 = 0.111690794839005; template<> Real GaussIntegrationTypeData<_git_triangle, 6>::quad_positions[] = {tri_6_b, tri_6_b, 1. - 2. * tri_6_b, tri_6_b, tri_6_b, 1. - 2. * tri_6_b, tri_6_a, 1. - 2. * tri_6_a, tri_6_a, tri_6_a, 1. - 2. * tri_6_a, tri_6_a}; template<> Real GaussIntegrationTypeData<_git_triangle, 6>::quad_weights[] = {tri_6_w1, tri_6_w1, tri_6_w1, tri_6_w2, tri_6_w2, tri_6_w2}; +/* -------------------------------------------------------------------------- */ +static const Real tri_7_a = (6. + std::sqrt(15.)) / 21.; +static const Real tri_7_b = (6. - std::sqrt(15.)) / 21.; +static const Real tri_7_w1 = (155. + std::sqrt(15.))/2400.; +static const Real tri_7_w2 = (155. - std::sqrt(15.))/2400.; +template<> Real GaussIntegrationTypeData<_git_triangle, 7>::quad_positions[] = { 1./3., 1./3., + tri_7_a, tri_7_a, + 1. - 2.*tri_7_a, tri_7_a, + tri_7_a, 1. - 2.*tri_7_a, + tri_7_b, tri_7_b, + 1. - 2.*tri_7_b, tri_7_b, + tri_7_b, 1. - 2.*tri_7_b}; +template<> Real GaussIntegrationTypeData<_git_triangle, 7>::quad_weights[] = {9./80., + tri_7_w1, tri_7_w1, tri_7_w1, + tri_7_w2, tri_7_w2, tri_7_w2}; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* Tetrahedrons */ /* -------------------------------------------------------------------------- */ template<> Real GaussIntegrationTypeData<_git_tetrahedron, 1>::quad_positions[] = {1./4., 1./4., 1./4.}; template<> Real GaussIntegrationTypeData<_git_tetrahedron, 1>::quad_weights[] = {1./6.}; /* -------------------------------------------------------------------------- */ static const Real tet_4_a = (5. - std::sqrt(5.))/20.; static const Real tet_4_b = (5. + 3.*std::sqrt(5.))/20.; template<> Real GaussIntegrationTypeData<_git_tetrahedron, 4>::quad_positions[] = {tet_4_a, tet_4_a, tet_4_a, tet_4_b, tet_4_a, tet_4_a, tet_4_a, tet_4_b, tet_4_a, tet_4_a, tet_4_a, tet_4_b}; template<> Real GaussIntegrationTypeData<_git_tetrahedron, 4>::quad_weights[] = {1./24., 1./24., 1./24., 1./24.}; /* -------------------------------------------------------------------------- */ template<> Real GaussIntegrationTypeData<_git_tetrahedron, 5>::quad_positions[] = {1./4., 1./4., 1./4., 1./6., 1./6., 1./6., 1./6., 1./6., 1./2., 1./6., 1./2., 1./6., 1./2., 1./6., 1./6.,}; template<> Real GaussIntegrationTypeData<_git_tetrahedron, 5>::quad_weights[] = {-2./15., 3./40., 3./40., 3./40., 3./40.}; /* -------------------------------------------------------------------------- */ static const Real tet_15_a = (7. + std::sqrt(15.))/34.; static const Real tet_15_b = (13. - 3. * std::sqrt(15.))/34.; static const Real tet_15_c = (7. - std::sqrt(15.))/34.; static const Real tet_15_d = (13. + 3. * std::sqrt(15.))/34.; static const Real tet_15_e = (5. - std::sqrt(15.))/20; static const Real tet_15_f = (5. + std::sqrt(15.))/20; static const Real tet_15_w1 = (2665. - 14. * std::sqrt(15.))/226800.; static const Real tet_15_w2 = (2665. + 14. * std::sqrt(15.))/226800.; static const Real tet_15_w3 = 5./567.; template<> Real GaussIntegrationTypeData<_git_tetrahedron, 15>::quad_positions[] = {1./4., 1./4., 1./4., tet_15_a, tet_15_a, tet_15_a, tet_15_a, tet_15_a, tet_15_b, tet_15_a, tet_15_b, tet_15_a, tet_15_b, tet_15_a, tet_15_a, tet_15_c, tet_15_c, tet_15_c, tet_15_c, tet_15_c, tet_15_d, tet_15_c, tet_15_d, tet_15_c, tet_15_d, tet_15_c, tet_15_c, tet_15_e, tet_15_e, tet_15_f, tet_15_e, tet_15_f, tet_15_e, tet_15_f, tet_15_e, tet_15_e, tet_15_e, tet_15_f, tet_15_f, tet_15_e, tet_15_f, tet_15_f, tet_15_f, tet_15_e, tet_15_f, tet_15_f, tet_15_f, tet_15_e}; template<> Real GaussIntegrationTypeData<_git_tetrahedron, 15>::quad_weights[] = {8./405., tet_15_w1, tet_15_w1, tet_15_w1, tet_15_w1, tet_15_w2, tet_15_w2, tet_15_w2, tet_15_w2, tet_15_w3, tet_15_w3, tet_15_w3, tet_15_w3, tet_15_w3, tet_15_w3}; -/* -------------------------------------------------------------------------- */ -/* Pentahedrons */ -/* -------------------------------------------------------------------------- */ -/// \todo write the cross product of _git_segment and _git_triangle -template<> Real GaussIntegrationTypeData<_git_pentahedron, 6>::quad_positions[] = {-1./sqrt(3.), 0.5, 0.5, - -1./sqrt(3.), 0. , 0.5, - -1./sqrt(3.), 0.5, 0., - 1./sqrt(3.), 0.5, 0.5, - 1./sqrt(3.), 0. , 0.5, - 1./sqrt(3.), 0.5 ,0.}; -template<> Real GaussIntegrationTypeData<_git_pentahedron, 6>::quad_weights[] = {1./6., 1./6., 1./6., - 1./6., 1./6., 1./6.}; -/* -------------------------------------------------------------------------- */ -template<> Real GaussIntegrationTypeData<_git_pentahedron, 8>::quad_positions[] = {-sqrt(3.)/3., 1./3., 1./3., - -sqrt(3.)/3., 0.6, 0.2, - -sqrt(3.)/3., 0.2, 0.6, - -sqrt(3.)/3., 0.2, 0.2, - sqrt(3.)/3., 1./3., 1./3., - sqrt(3.)/3., 0.6, 0.2, - sqrt(3.)/3., 0.2, 0.6, - sqrt(3.)/3., 0.2, 0.2}; -template<> Real GaussIntegrationTypeData<_git_pentahedron, 8>::quad_weights[] = {-27./96., 25./96., 25./96., 25./96., - -27./96., 25./96., 25./96., 25./96.}; -/* -------------------------------------------------------------------------- */ -static const Real pent_21_x = std::sqrt(3./5.); -static const Real pent_21_a = (6. + std::sqrt(15.)) / 21.; -static const Real pent_21_b = (6. - std::sqrt(15.)) / 21.; -static const Real pent_21_w1_1 = 5./9.; -static const Real pent_21_w2_1 = 8./9.; -static const Real pent_21_w1_2 = (155. + std::sqrt(15.))/2400.; -static const Real pent_21_w2_2 = (155. - std::sqrt(15.))/2400.; -template<> Real GaussIntegrationTypeData<_git_pentahedron, 21>::quad_positions[] = {- pent_21_x, 1./3., 1./3., - - pent_21_x, pent_21_a, pent_21_a, - - pent_21_x, 1. - 2.*pent_21_a, pent_21_a, - - pent_21_x, pent_21_a, 1. - 2.*pent_21_a, - - pent_21_x, pent_21_b, pent_21_b, - - pent_21_x, 1. - 2.*pent_21_b, pent_21_b, - - pent_21_x, pent_21_b, 1. - 2.*pent_21_b, - 0., 1./3., 1./3., - 0., pent_21_a, pent_21_a, - 0., 1. - 2.*pent_21_a, pent_21_a, - 0., pent_21_a, 1. - 2.*pent_21_a, - 0., pent_21_b, pent_21_b, - 0., 1. - 2.*pent_21_b, pent_21_b, - 0., pent_21_b, 1. - 2.*pent_21_b, - pent_21_x, 1./3., 1./3., - pent_21_x, pent_21_a, pent_21_a, - pent_21_x, 1. - 2.*pent_21_a, pent_21_a, - pent_21_x, pent_21_a, 1. - 2.*pent_21_a, - pent_21_x, pent_21_b, pent_21_b, - pent_21_x, 1. - 2.*pent_21_b, pent_21_b, - pent_21_x, pent_21_b, 1. - 2.*pent_21_b}; -template<> Real GaussIntegrationTypeData<_git_pentahedron, 21>::quad_weights[] = {1./16., - pent_21_w1_1*pent_21_w1_2, pent_21_w1_1*pent_21_w1_2, pent_21_w1_1*pent_21_w1_2, - pent_21_w1_1*pent_21_w2_2, pent_21_w1_1*pent_21_w2_2, pent_21_w1_1*pent_21_w2_2, - 1./10., - pent_21_w1_2*pent_21_w1_2, pent_21_w1_2*pent_21_w1_2, pent_21_w1_2*pent_21_w1_2, - pent_21_w1_2*pent_21_w2_2, pent_21_w1_2*pent_21_w2_2, pent_21_w1_2*pent_21_w2_2, - 1./16., - pent_21_w1_1*pent_21_w1_2, pent_21_w1_1*pent_21_w1_2, pent_21_w1_1*pent_21_w1_2, - pent_21_w1_1*pent_21_w2_2, pent_21_w1_1*pent_21_w2_2, pent_21_w1_1*pent_21_w2_2}; __END_AKANTU__ diff --git a/src/fe_engine/gauss_integration_tmpl.hh b/src/fe_engine/gauss_integration_tmpl.hh new file mode 100644 index 000000000..2e7c36460 --- /dev/null +++ b/src/fe_engine/gauss_integration_tmpl.hh @@ -0,0 +1,289 @@ +/** + * @file gauss_integration_tmpl.hh + * + * @author Nicolas Richart + * + * @date Tue May 10 10:48:21 2016 + * + * @brief implementation of the gauss integration helpers + * + * @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_GAUSS_INTEGRATION_TMPL_HH__ +#define __AKANTU_GAUSS_INTEGRATION_TMPL_HH__ + +__BEGIN_AKANTU__ + +/* -------------------------------------------------------------------------- */ +/* GaussIntegrationElement */ +/* -------------------------------------------------------------------------- */ +namespace _aka_gauss_helpers { +template struct GaussIntegrationNbPoints { + static const UInt nb_points = -1; +}; + +template struct GaussIntegrationNbPoints<_git_point, n> { + static const UInt nb_points = 1; +}; + +template struct GaussIntegrationNbPoints<_git_segment, n> { + static const UInt nb_points = (n + 1) / 2 + ((n + 1) % 2 ? 1 : 0); +}; + +#define DECLARE_GAUSS_NB_POINTS(type, order, points) \ + template <> struct GaussIntegrationNbPoints { \ + static const UInt nb_points = points; \ + } + +#define DECLARE_GAUSS_NB_POINTS_PENT(type, order, xo, yo) \ + template <> struct GaussIntegrationNbPoints { \ + static const UInt x_order = xo; \ + static const UInt yz_order = yo; \ + static const UInt nb_points = 1; \ + } + +DECLARE_GAUSS_NB_POINTS(_git_triangle, 1, 1); +DECLARE_GAUSS_NB_POINTS(_git_triangle, 2, 3); +DECLARE_GAUSS_NB_POINTS(_git_triangle, 3, 4); +DECLARE_GAUSS_NB_POINTS(_git_triangle, 4, 6); +DECLARE_GAUSS_NB_POINTS(_git_triangle, 5, 7); +DECLARE_GAUSS_NB_POINTS(_git_tetrahedron, 1, 1); +DECLARE_GAUSS_NB_POINTS(_git_tetrahedron, 2, 4); +DECLARE_GAUSS_NB_POINTS(_git_tetrahedron, 3, 5); +DECLARE_GAUSS_NB_POINTS(_git_tetrahedron, 4, 15); +DECLARE_GAUSS_NB_POINTS(_git_tetrahedron, 5, 15); +DECLARE_GAUSS_NB_POINTS_PENT(_git_pentahedron, 1, 3, + 2); // order 3 in x, order 2 in y and z +DECLARE_GAUSS_NB_POINTS_PENT(_git_pentahedron, 2, 3, + 2); // order 3 in x, order 2 in y and z +DECLARE_GAUSS_NB_POINTS_PENT(_git_pentahedron, 3, 3, + 3); // order 3 in x, order 3 in y and z +DECLARE_GAUSS_NB_POINTS_PENT(_git_pentahedron, 4, 5, + 5); // order 5 in x, order 5 in y and z +DECLARE_GAUSS_NB_POINTS_PENT(_git_pentahedron, 5, 5, + 5); // order 5 in x, order 5 in y and z + +template +struct GaussIntegrationNbPointsHelper { + static const UInt pnp = GaussIntegrationNbPoints::nb_points; + static const UInt order = n; + static const UInt nb_points = pnp; +}; + +template +struct GaussIntegrationNbPointsHelper { + static const UInt nb_points = 0; +}; + +template +struct GaussIntegrationTypeMissingData { + typedef void value; +}; + +template +struct GaussIntegrationTypeMissingData { + // If you fall here it means you are missing GaussIntegrationTypeData for + // the GaussIntegrationType 'type' with polynomials of order 'n' +}; + +template +struct GaussIntegrationTypeMissingData { + // If you fall here it means you are missing GaussIntegrationTypeData for + // the GaussIntegrationType 'type' with polynomials of order 'n' +}; + +/* ------------------------------------------------------------------------ */ +/* Generic helper */ +/* ------------------------------------------------------------------------ */ +template +struct GaussIntegrationTypeDataHelper { + typedef GaussIntegrationNbPoints git_np; + typedef GaussIntegrationTypeData git_data; + typedef typename GaussIntegrationTypeMissingData< + type, n, git_np::nb_points>::value git_missing; + + static UInt getNbQuadraturePoints() { return git_np::nb_points; } + + static const Matrix getQuadraturePoints() { + return Matrix(git_data::quad_positions, dimension, git_np::nb_points); + } + + static const Vector getWeights() { + return Vector(git_data::quad_weights, git_np::nb_points); + } +}; + +/* ------------------------------------------------------------------------ */ +/* helper for _segment _quadrangle _hexahedron */ +/* ------------------------------------------------------------------------ */ +template +struct GaussIntegrationTypeDataHelper<_git_segment, dimension, dp> { + typedef GaussIntegrationNbPoints<_git_segment, dp> git_np; + typedef GaussIntegrationTypeData<_git_segment, git_np::nb_points> git_data; + typedef typename GaussIntegrationTypeMissingData< + _git_segment, dp, git_np::nb_points>::value git_missing; + + static UInt getNbQuadraturePoints() { + return Math::pow(git_np::nb_points); + } + + static const Matrix getQuadraturePoints() { + UInt tot_nquad = getNbQuadraturePoints(); + UInt nquad = git_np::nb_points; + + Matrix quads(dimension, tot_nquad); + Vector pos(git_data::quad_positions, nquad); + + UInt offset = 1; + for (UInt d = 0; d < dimension; ++d) { + for (UInt n = 0, q = 0; n < tot_nquad; ++n, q += offset) { + UInt rq = q % tot_nquad + q / tot_nquad; + quads(d, rq) = pos(n % nquad); + } + offset *= nquad; + } + return quads; + } + + static const Vector getWeights() { + UInt tot_nquad = getNbQuadraturePoints(); + UInt nquad = git_np::nb_points; + + Vector quads_weights(tot_nquad, 1.); + Vector weights(git_data::quad_weights, nquad); + + UInt offset = 1; + for (UInt d = 0; d < dimension; ++d) { + for (UInt n = 0, q = 0; n < tot_nquad; ++n, q += offset) { + UInt rq = q % tot_nquad + q / tot_nquad; + quads_weights(rq) *= weights(n % nquad); + } + offset *= nquad; + } + return quads_weights; + } +}; + +/* ------------------------------------------------------------------------ */ +/* helper for _pentahedron */ +/* ------------------------------------------------------------------------ */ +template +struct GaussIntegrationTypeDataHelper<_git_pentahedron, dimension, dp> { + typedef GaussIntegrationNbPoints<_git_pentahedron, dp> git_info; + typedef GaussIntegrationNbPoints<_git_segment, git_info::x_order> git_np_seg; + typedef GaussIntegrationNbPoints<_git_triangle, git_info::yz_order> git_np_tri; + + typedef GaussIntegrationTypeData<_git_segment, git_np_seg::nb_points> + git_data_seg; + typedef GaussIntegrationTypeData<_git_triangle, git_np_tri::nb_points> + git_data_tri; + typedef typename GaussIntegrationTypeMissingData<_git_segment, git_info::x_order, + git_np_seg::nb_points>::value git_missing_seg; + typedef typename GaussIntegrationTypeMissingData<_git_triangle, git_info::yz_order, + git_np_tri::nb_points>::value git_missing_tri; + + static UInt getNbQuadraturePoints() { + return git_np_seg::nb_points * git_np_tri::nb_points; + } + + static const Matrix getQuadraturePoints() { + UInt tot_nquad = getNbQuadraturePoints(); + UInt nquad_seg = git_np_seg::nb_points; + UInt nquad_tri = git_np_tri::nb_points; + + Matrix quads(dimension, tot_nquad); + Vector pos_seg(git_data_seg::quad_positions, nquad_seg); + Matrix pos_tri(git_data_tri::quad_positions, 2, nquad_tri); + + for (UInt ns = 0, q = 0; ns < nquad_seg; ++ns) { + for (UInt nt = 0; nt < nquad_tri; ++nt, ++q) { + quads(q, _x) = pos_seg(ns); + quads(q, _y) = pos_tri(nt, _x); + quads(q, _z) = pos_tri(nt, _y); + } + } + return quads; + } + + static const Vector getWeights() { + UInt tot_nquad = getNbQuadraturePoints(); + UInt nquad_seg = git_np_seg::nb_points; + UInt nquad_tri = git_np_tri::nb_points; + + Vector quads_weights(tot_nquad); + + Vector weight_seg(git_data_seg::quad_weights, nquad_seg); + Vector weight_tri(git_data_tri::quad_weights, nquad_tri); + + for (UInt ns = 0, q = 0; ns < nquad_seg; ++ns) { + for (UInt nt = 0; nt < nquad_tri; ++nt, ++q) { + quads_weights(q) = weight_seg(ns) * weight_tri(nt); + } + } + return quads_weights; + } +}; +} + +template +const Matrix +GaussIntegrationElement::getQuadraturePoints() { + const InterpolationType itp_type = + ElementClassProperty::interpolation_type; + typedef InterpolationPorperty interpolation_property; + typedef _aka_gauss_helpers::GaussIntegrationTypeDataHelper< + ElementClassProperty::gauss_integration_type, + interpolation_property::natural_space_dimension, n> data_helper; + Matrix tmp(data_helper::getQuadraturePoints()); + return tmp; +} + +/* -------------------------------------------------------------------------- */ +template +const Vector GaussIntegrationElement::getWeights() { + const InterpolationType itp_type = + ElementClassProperty::interpolation_type; + typedef InterpolationPorperty interpolation_property; + typedef _aka_gauss_helpers::GaussIntegrationTypeDataHelper< + ElementClassProperty::gauss_integration_type, + interpolation_property::natural_space_dimension, n> data_helper; + Vector tmp(data_helper::getWeights()); + return tmp; +} + +/* -------------------------------------------------------------------------- */ +template +UInt GaussIntegrationElement::getNbQuadraturePoints() { + const InterpolationType itp_type = + ElementClassProperty::interpolation_type; + typedef InterpolationPorperty interpolation_property; + typedef _aka_gauss_helpers::GaussIntegrationTypeDataHelper< + ElementClassProperty::gauss_integration_type, + interpolation_property::natural_space_dimension, n> data_helper; + return data_helper::getNbQuadraturePoints(); +} + +__END_AKANTU__ + +#endif /* __AKANTU_GAUSS_INTEGRATION_TMPL_HH__ */