diff --git a/packages/core.cmake b/packages/core.cmake index c0a9178ea..24b09898a 100644 --- a/packages/core.cmake +++ b/packages/core.cmake @@ -1,529 +1,534 @@ #=============================================================================== # @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" FEATURES_PUBLIC cxx_strong_enums cxx_defaulted_functions cxx_deleted_functions cxx_auto_type cxx_decltype_auto FEATURES_PRIVATE cxx_lambdas cxx_nullptr cxx_range_for cxx_delegating_constructors DEPENDS INTERFACE akantu_iterators Boost) package_declare_sources(core common/aka_array.cc common/aka_array.hh common/aka_array_tmpl.hh common/aka_array_printer.hh common/aka_bbox.hh common/aka_blas_lapack.hh common/aka_circular_array.hh common/aka_circular_array_inline_impl.hh common/aka_common.cc common/aka_common.hh common/aka_common_inline_impl.hh common/aka_csr.hh common/aka_element_classes_info_inline_impl.hh common/aka_enum_macros.hh common/aka_error.cc common/aka_error.hh common/aka_event_handler_manager.hh common/aka_extern.cc common/aka_factory.hh 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.hh common/aka_named_argument.hh 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.hh 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_tmpl.hh common/aka_voigthelper.cc common/aka_warning.hh common/aka_warning_restore.hh fe_engine/element_class.hh + fe_engine/element_class_helper.hh fe_engine/element_class_tmpl.hh fe_engine/element_classes/element_class_hexahedron_8_inline_impl.hh fe_engine/element_classes/element_class_hexahedron_20_inline_impl.hh fe_engine/element_classes/element_class_pentahedron_6_inline_impl.hh fe_engine/element_classes/element_class_pentahedron_15_inline_impl.hh fe_engine/element_classes/element_class_point_1_inline_impl.hh fe_engine/element_classes/element_class_quadrangle_4_inline_impl.hh fe_engine/element_classes/element_class_quadrangle_8_inline_impl.hh fe_engine/element_classes/element_class_segment_2_inline_impl.hh fe_engine/element_classes/element_class_segment_3_inline_impl.hh fe_engine/element_classes/element_class_tetrahedron_10_inline_impl.hh fe_engine/element_classes/element_class_tetrahedron_4_inline_impl.hh fe_engine/element_classes/element_class_triangle_3_inline_impl.hh fe_engine/element_classes/element_class_triangle_6_inline_impl.hh fe_engine/element_type_conversion.hh fe_engine/fe_engine.cc fe_engine/fe_engine.hh fe_engine/fe_engine_inline_impl.hh fe_engine/fe_engine_template.hh fe_engine/fe_engine_template_tmpl_field.hh fe_engine/fe_engine_template_tmpl.hh fe_engine/geometrical_element_property.hh fe_engine/geometrical_element_property.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.hh fe_engine/interpolation_element_tmpl.hh fe_engine/integration_point.hh fe_engine/shape_functions.hh fe_engine/shape_functions.cc fe_engine/shape_functions_inline_impl.hh fe_engine/shape_lagrange_base.cc fe_engine/shape_lagrange_base.hh fe_engine/shape_lagrange_base_inline_impl.hh fe_engine/shape_lagrange.hh fe_engine/shape_lagrange_inline_impl.hh 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_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/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 io/parser/parameter_registry.cc io/parser/parameter_registry.hh io/parser/parameter_registry_tmpl.hh mesh/element_group.cc mesh/element_group.hh mesh/element_group_inline_impl.hh mesh/element_type_map.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.hh mesh/mesh.cc mesh/mesh.hh mesh/mesh_periodic.cc mesh/mesh_accessor.hh mesh/mesh_events.hh mesh/mesh_filter.hh mesh/mesh_global_data_updater.hh mesh/mesh_data.cc mesh/mesh_data.hh mesh/mesh_data_tmpl.hh mesh/mesh_inline_impl.hh mesh/node_group.cc mesh/node_group.hh mesh/node_group_inline_impl.hh mesh/mesh_iterators.hh 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_distribution.cc mesh_utils/mesh_utils_distribution.hh mesh_utils/mesh_utils.hh mesh_utils/mesh_utils_inline_impl.hh mesh_utils/global_ids_updater.hh mesh_utils/global_ids_updater.cc mesh_utils/global_ids_updater_inline_impl.hh model/common/boundary_condition/boundary_condition.hh model/common/boundary_condition/boundary_condition_functor.hh model/common/boundary_condition/boundary_condition_functor_inline_impl.hh model/common/boundary_condition/boundary_condition_tmpl.hh model/common/non_local_toolbox/neighborhood_base.hh model/common/non_local_toolbox/neighborhood_base.cc model/common/non_local_toolbox/neighborhood_base_inline_impl.hh model/common/non_local_toolbox/neighborhoods_criterion_evaluation/neighborhood_max_criterion.hh model/common/non_local_toolbox/neighborhoods_criterion_evaluation/neighborhood_max_criterion.cc model/common/non_local_toolbox/neighborhoods_criterion_evaluation/neighborhood_max_criterion_inline_impl.hh model/common/non_local_toolbox/non_local_manager.hh model/common/non_local_toolbox/non_local_manager.cc model/common/non_local_toolbox/non_local_manager_inline_impl.hh model/common/non_local_toolbox/non_local_manager_callback.hh model/common/non_local_toolbox/non_local_neighborhood_base.hh model/common/non_local_toolbox/non_local_neighborhood_base.cc model/common/non_local_toolbox/non_local_neighborhood.hh model/common/non_local_toolbox/non_local_neighborhood_tmpl.hh model/common/non_local_toolbox/non_local_neighborhood_inline_impl.hh model/common/non_local_toolbox/base_weight_function.hh model/common/non_local_toolbox/base_weight_function_inline_impl.hh model/common/model_solver.cc model/common/model_solver.hh model/common/solver_callback.hh model/common/solver_callback.cc model/common/dof_manager/dof_manager.cc model/common/dof_manager/dof_manager.hh model/common/dof_manager/dof_manager_default.cc model/common/dof_manager/dof_manager_default.hh model/common/dof_manager/dof_manager_default_inline_impl.hh model/common/dof_manager/dof_manager_inline_impl.hh model/common/non_linear_solver/non_linear_solver.cc model/common/non_linear_solver/non_linear_solver.hh model/common/non_linear_solver/non_linear_solver_default.hh model/common/non_linear_solver/non_linear_solver_lumped.cc model/common/non_linear_solver/non_linear_solver_lumped.hh model/common/time_step_solvers/time_step_solver.hh model/common/time_step_solvers/time_step_solver.cc model/common/time_step_solvers/time_step_solver_default.cc model/common/time_step_solvers/time_step_solver_default.hh model/common/time_step_solvers/time_step_solver_default_explicit.hh model/common/integration_scheme/generalized_trapezoidal.cc model/common/integration_scheme/generalized_trapezoidal.hh model/common/integration_scheme/integration_scheme.cc model/common/integration_scheme/integration_scheme.hh model/common/integration_scheme/integration_scheme_1st_order.cc model/common/integration_scheme/integration_scheme_1st_order.hh model/common/integration_scheme/integration_scheme_2nd_order.cc model/common/integration_scheme/integration_scheme_2nd_order.hh model/common/integration_scheme/newmark-beta.cc model/common/integration_scheme/newmark-beta.hh model/common/integration_scheme/pseudo_time.cc model/common/integration_scheme/pseudo_time.hh model/model.cc model/model.hh model/model_inline_impl.hh model/model_options.hh solver/solver_vector.hh solver/solver_vector_default.cc solver/solver_vector_default.hh solver/solver_vector_default_tmpl.hh solver/solver_vector_distributed.cc solver/solver_vector_distributed.hh solver/sparse_matrix.cc solver/sparse_matrix.hh solver/sparse_matrix_aij.cc solver/sparse_matrix_aij.hh solver/sparse_matrix_aij_inline_impl.hh solver/sparse_matrix_inline_impl.hh solver/sparse_solver.cc solver/sparse_solver.hh solver/sparse_solver_inline_impl.hh solver/terms_to_assemble.hh synchronizer/communication_buffer_inline_impl.hh synchronizer/communication_descriptor.hh synchronizer/communication_descriptor_tmpl.hh synchronizer/communication_request.hh synchronizer/communication_tag.hh synchronizer/communications.hh synchronizer/communications_tmpl.hh synchronizer/communicator.cc synchronizer/communicator.hh synchronizer/communicator_dummy_inline_impl.hh synchronizer/communicator_event_handler.hh synchronizer/communicator_inline_impl.hh synchronizer/data_accessor.cc synchronizer/data_accessor.hh synchronizer/dof_synchronizer.cc synchronizer/dof_synchronizer.hh synchronizer/dof_synchronizer_inline_impl.hh synchronizer/element_info_per_processor.cc synchronizer/element_info_per_processor.hh synchronizer/element_info_per_processor_tmpl.hh synchronizer/element_synchronizer.cc synchronizer/element_synchronizer.hh synchronizer/facet_synchronizer.cc synchronizer/facet_synchronizer.hh synchronizer/facet_synchronizer_inline_impl.hh synchronizer/grid_synchronizer.cc synchronizer/grid_synchronizer.hh synchronizer/grid_synchronizer_tmpl.hh synchronizer/master_element_info_per_processor.cc synchronizer/node_info_per_processor.cc synchronizer/node_info_per_processor.hh synchronizer/node_synchronizer.cc synchronizer/node_synchronizer.hh synchronizer/node_synchronizer_inline_impl.hh synchronizer/periodic_node_synchronizer.cc synchronizer/periodic_node_synchronizer.hh synchronizer/slave_element_info_per_processor.cc synchronizer/synchronizer.cc synchronizer/synchronizer.hh synchronizer/synchronizer_impl.hh synchronizer/synchronizer_impl_tmpl.hh synchronizer/synchronizer_registry.cc synchronizer/synchronizer_registry.hh synchronizer/synchronizer_tmpl.hh synchronizer/communication_buffer.hh ) set(AKANTU_SPIRIT_SOURCES io/mesh_io/mesh_io_abaqus.cc io/parser/parser_real.cc io/parser/parser_random.cc io/parser/parser_types.cc io/parser/parser_input_files.cc PARENT_SCOPE ) 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_N + compute_dnds + compute_d2nds2 + compute_jmat get_shapes_derivatives lagrange_base assemble_fields ) package_declare_documentation_files(core manual.sty manual.cls manual.tex manual-macros.sty manual-titlepages.tex manual-authors.tex manual-changelog.tex manual-introduction.tex manual-gettingstarted.tex manual-io.tex manual-feengine.tex manual-elements.tex manual-appendix-elements.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/hot-point-1.png figures/hot-point-2.png figures/insertion.pdf figures/interpolate.pdf figures/interpolate.svg figures/vectors.pdf figures/vectors.svg 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) package_declare_extra_files_to_package(core SOURCES common/aka_element_classes_info.hh.in common/aka_config.hh.in ) if(CMAKE_CXX_COMPILER_ID STREQUAL "Clang" AND (NOT CMAKE_CXX_COMPILER_VERSION VERSION_LESS 3.9)) package_set_compile_flags(core CXX "-Wno-undefined-var-template") endif() if(DEFINED AKANTU_CXX11_FLAGS) package_declare(core_cxx11 NOT_OPTIONAL DESCRIPTION "C++ 11 additions for Akantu core" COMPILE_FLAGS CXX "${AKANTU_CXX11_FLAGS}") if(CMAKE_CXX_COMPILER_ID STREQUAL "GNU") if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS "4.6") set(AKANTU_CORE_CXX11 OFF CACHE BOOL "C++ 11 additions for Akantu core - not supported by the selected compiler" FORCE) endif() endif() package_declare_documentation(core_cxx11 "This option activates some features of the C++11 standard. This is usable with GCC>=4.7 or Intel>=13.") else() if(CMAKE_VERSION VERSION_LESS 3.1) message(FATAL_ERROR "Since version 3.0 Akantu requires at least c++11 capable compiler") endif() endif() diff --git a/src/fe_engine/element_class.hh b/src/fe_engine/element_class.hh index 51aabf127..c8b742284 100644 --- a/src/fe_engine/element_class.hh +++ b/src/fe_engine/element_class.hh @@ -1,432 +1,432 @@ /** * @file element_class.hh * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Feb 20 2018 * * @brief Declaration of the ElementClass main class and the * Integration and Interpolation elements * * * Copyright (©) 2010-2018 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_types.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_ELEMENT_CLASS_HH_ #define AKANTU_ELEMENT_CLASS_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ /// default element class structure template struct ElementClassProperty { static const GeometricalType geometrical_type{_gt_not_defined}; static const InterpolationType interpolation_type{_itp_not_defined}; static const ElementKind element_kind{_ek_regular}; static const UInt spatial_dimension{0}; static const GaussIntegrationType gauss_integration_type{_git_not_defined}; static const UInt polynomial_degree{0}; }; #if !defined(DOXYGEN) /// Macro to generate the element class structures for different element types #define AKANTU_DEFINE_ELEMENT_CLASS_PROPERTY(elem_type, geom_type, \ interp_type, elem_kind, sp, \ gauss_int_type, min_int_order) \ template <> struct ElementClassProperty { \ static const GeometricalType geometrical_type{geom_type}; \ static const InterpolationType interpolation_type{interp_type}; \ static const ElementKind element_kind{elem_kind}; \ static const UInt spatial_dimension{sp}; \ static const GaussIntegrationType gauss_integration_type{gauss_int_type}; \ static const UInt polynomial_degree{min_int_order}; \ } #else #define AKANTU_DEFINE_ELEMENT_CLASS_PROPERTY(elem_type, geom_type, \ interp_type, elem_kind, sp, \ gauss_int_type, min_int_order) #endif /* -------------------------------------------------------------------------- */ /* Geometry */ /* -------------------------------------------------------------------------- */ /// Default GeometricalShape structure template struct GeometricalShape { static const GeometricalShapeType shape{_gst_point}; }; /// Templated GeometricalShape with function contains template struct GeometricalShapeContains { /// Check if the point (vector in 2 and 3D) at natural coordinate coor template static inline bool contains(const vector_type & coord); }; #if !defined(DOXYGEN) /// Macro to generate the GeometricalShape structures for different geometrical /// types #define AKANTU_DEFINE_SHAPE(geom_type, geom_shape) \ template <> struct GeometricalShape { \ static const GeometricalShapeType shape{geom_shape}; \ } AKANTU_DEFINE_SHAPE(_gt_hexahedron_20, _gst_square); AKANTU_DEFINE_SHAPE(_gt_hexahedron_8, _gst_square); AKANTU_DEFINE_SHAPE(_gt_pentahedron_15, _gst_prism); AKANTU_DEFINE_SHAPE(_gt_pentahedron_6, _gst_prism); AKANTU_DEFINE_SHAPE(_gt_point, _gst_point); AKANTU_DEFINE_SHAPE(_gt_quadrangle_4, _gst_square); AKANTU_DEFINE_SHAPE(_gt_quadrangle_8, _gst_square); AKANTU_DEFINE_SHAPE(_gt_segment_2, _gst_square); AKANTU_DEFINE_SHAPE(_gt_segment_3, _gst_square); AKANTU_DEFINE_SHAPE(_gt_tetrahedron_10, _gst_triangle); AKANTU_DEFINE_SHAPE(_gt_tetrahedron_4, _gst_triangle); AKANTU_DEFINE_SHAPE(_gt_triangle_3, _gst_triangle); AKANTU_DEFINE_SHAPE(_gt_triangle_6, _gst_triangle); #endif /* -------------------------------------------------------------------------- */ template struct GeometricalElementProperty {}; template struct ElementClassExtraGeometryProperties {}; /* -------------------------------------------------------------------------- */ /// Templated GeometricalElement with function getInradius template ::shape> class GeometricalElement { using geometrical_property = GeometricalElementProperty; public: /// compute the in-radius: \todo should be renamed for characteristic length static inline Real getInradius(const Matrix & /*coord*/) { AKANTU_TO_IMPLEMENT(); } /// true if the natural coordinates are in the element template static inline bool contains(const vector_type & coord); public: static AKANTU_GET_MACRO_NOT_CONST(SpatialDimension, geometrical_property::spatial_dimension, UInt); static AKANTU_GET_MACRO_NOT_CONST(NbNodesPerElement, geometrical_property::nb_nodes_per_element, UInt); static inline constexpr auto getNbFacetTypes() { return geometrical_property::nb_facet_types; }; static inline UInt getNbFacetsPerElement(UInt t); static inline UInt getNbFacetsPerElement(); static inline constexpr auto getFacetLocalConnectivityPerElement(UInt t = 0); }; /* -------------------------------------------------------------------------- */ /* Interpolation */ /* -------------------------------------------------------------------------- */ /// default InterpolationProperty structure template struct InterpolationProperty {}; #if !defined(DOXYGEN) /// Macro to generate the InterpolationProperty structures for different /// interpolation types #define AKANTU_DEFINE_INTERPOLATION_TYPE_PROPERTY(itp_type, itp_kind, \ nb_nodes, ndim) \ template <> struct InterpolationProperty { \ static constexpr InterpolationKind kind{itp_kind}; \ static constexpr UInt nb_nodes_per_element{nb_nodes}; \ static constexpr UInt natural_space_dimension{ndim}; \ } #else #define AKANTU_DEFINE_INTERPOLATION_TYPE_PROPERTY(itp_type, itp_kind, \ nb_nodes, ndim) #endif /* -------------------------------------------------------------------------- */ /// Generic (templated by the enum InterpolationType which specifies the order /// and the dimension of the interpolation) class handling the elemental /// interpolation template ::kind> class InterpolationElement { public: using interpolation_property = InterpolationProperty; /// compute the shape values for a given set of points in natural coordinates static inline void computeShapes(const Matrix & natural_coord, Matrix & N); /// compute the shape values for a given point in natural coordinates template static inline void computeShapes(const vector_type & /*unused*/, vector_type & /*unused*/) { AKANTU_TO_IMPLEMENT(); } /** * compute @f$ B_{ij} = \frac{\partial N_j}{\partial S_i} @f$ the variation of * shape functions along with variation of natural coordinates on a given set * of points in natural coordinates */ static inline void computeDNDS(const Matrix & natural_coord, Tensor3 & dnds); /** * compute @f$ B_{ij} = \frac{\partial N_j}{\partial S_i} @f$ the variation of * shape functions along with * variation of natural coordinates on a given point in natural * coordinates */ template static inline void computeDNDS(const vector_type & /*unused*/, matrix_type & /*unused*/) { AKANTU_TO_IMPLEMENT(); } /** * compute @f$ @f$ **/ - static inline void computeDN2DS2(const Matrix & natural_coord, - Tensor3 & dn2ds2); + static inline void computeD2NDS2(const Matrix & natural_coord, + Tensor3 & d2nds2); /** * compute @f$ B_{ij} = \frac{\partial N_j}{\partial S_i} @f$ the * second variation of * shape functions along with * variation of natural coordinates on a given point in natural * coordinates */ template - static inline void computeDN2DS2(const vector_type &, matrix_type &) { + static inline void computeD2NDS2(const vector_type &, matrix_type &) { AKANTU_TO_IMPLEMENT(); } /// compute jacobian (or integration variable change factor) for a given point /// in the case of spatial_dimension != natural_space_dimension static inline void computeSpecialJacobian(const Matrix & /*unused*/, Real & /*unused*/) { AKANTU_TO_IMPLEMENT(); } /// interpolate a field given (arbitrary) natural coordinates static inline void interpolateOnNaturalCoordinates(const Vector & natural_coords, const Matrix & nodal_values, Vector & interpolated); /// interpolate a field given the shape functions on the interpolation point static inline void interpolate(const Matrix & nodal_values, const Vector & shapes, Vector & interpolated); /// interpolate a field given the shape functions on the interpolations points static inline void interpolate(const Matrix & nodal_values, const Matrix & shapes, Matrix & interpolated); /// compute the gradient of a given field on the given natural coordinates static inline void gradientOnNaturalCoordinates(const Vector & natural_coords, const Matrix & f, Matrix & gradient); public: static AKANTU_GET_MACRO_NOT_CONST( ShapeSize, InterpolationProperty::nb_nodes_per_element, UInt); static AKANTU_GET_MACRO_NOT_CONST( ShapeDerivativesSize, (InterpolationProperty::nb_nodes_per_element * InterpolationProperty::natural_space_dimension), UInt); static AKANTU_GET_MACRO_NOT_CONST( NaturalSpaceDimension, InterpolationProperty::natural_space_dimension, UInt); static AKANTU_GET_MACRO_NOT_CONST( NbNodesPerInterpolationElement, InterpolationProperty::nb_nodes_per_element, UInt); }; /* -------------------------------------------------------------------------- */ /* Integration */ /* -------------------------------------------------------------------------- */ template struct GaussIntegrationTypeData { /// quadrature points in natural coordinates static Real quad_positions[]; /// weights for the Gauss integration static Real quad_weights[]; }; template ::polynomial_degree> class GaussIntegrationElement { public: static UInt getNbQuadraturePoints(); static Matrix getQuadraturePoints(); static Vector getWeights(); }; /* -------------------------------------------------------------------------- */ /* ElementClass */ /* -------------------------------------------------------------------------- */ template ::element_kind> class ElementClass : public GeometricalElement< ElementClassProperty::geometrical_type>, public InterpolationElement< ElementClassProperty::interpolation_type> { protected: using geometrical_element = GeometricalElement::geometrical_type>; using interpolation_element = InterpolationElement< ElementClassProperty::interpolation_type>; using element_property = ElementClassProperty; using interpolation_property = typename interpolation_element::interpolation_property; public: /** * compute @f$ J = \frac{\partial x_j}{\partial s_i} @f$ the variation of real * coordinates along with variation of natural coordinates on a given point in * natural coordinates */ static inline void computeJMat(const Matrix & dnds, const Matrix & node_coords, Matrix & J); /** * compute the Jacobian matrix by computing the variation of real coordinates * along with variation of natural coordinates on a given set of points in * natural coordinates */ static inline void computeJMat(const Tensor3 & dnds, const Matrix & node_coords, Tensor3 & J); /// compute the jacobians of a serie of natural coordinates static inline void computeJacobian(const Matrix & natural_coords, const Matrix & node_coords, Vector & jacobians); /// compute jacobian (or integration variable change factor) for a set of /// points static inline void computeJacobian(const Tensor3 & J, Vector & jacobians); /// compute jacobian (or integration variable change factor) for a given point static inline void computeJacobian(const Matrix & J, Real & jacobians); /// compute shape derivatives (input is dxds) for a set of points static inline void computeShapeDerivatives(const Tensor3 & J, const Tensor3 & dnds, Tensor3 & shape_deriv); /// compute shape derivatives (input is dxds) for a given point static inline void computeShapeDerivatives(const Matrix & J, const Matrix & dnds, Matrix & shape_deriv); /// compute the normal of a surface defined by the function f static inline void computeNormalsOnNaturalCoordinates(const Matrix & coord, Matrix & f, Matrix & normals); /// get natural coordinates from real coordinates static inline void inverseMap(const Vector & real_coords, const Matrix & node_coords, Vector & natural_coords, UInt max_iterations = 100, Real tolerance = 1e-10); /// get natural coordinates from real coordinates static inline void inverseMap(const Matrix & real_coords, const Matrix & node_coords, Matrix & natural_coords, UInt max_iterations = 100, Real tolerance = 1e-10); public: static AKANTU_GET_MACRO_NOT_CONST(Kind, element_kind, ElementKind); static constexpr AKANTU_GET_MACRO_NOT_CONST( SpatialDimension, ElementClassProperty::spatial_dimension, UInt); using element_class_extra_geom_property = ElementClassExtraGeometryProperties; static constexpr auto getP1ElementType() { return element_class_extra_geom_property::p1_type; } static constexpr auto getFacetType(UInt t = 0) { return element_class_extra_geom_property::facet_type[t]; } static constexpr auto getFacetTypes(); }; /* -------------------------------------------------------------------------- */ } // namespace akantu /* -------------------------------------------------------------------------- */ #include "geometrical_element_property.hh" #include "interpolation_element_tmpl.hh" /* -------------------------------------------------------------------------- */ #include "element_class_tmpl.hh" /* -------------------------------------------------------------------------- */ #include "element_class_hexahedron_8_inline_impl.hh" #include "element_class_pentahedron_6_inline_impl.hh" /* keep order */ #include "element_class_hexahedron_20_inline_impl.hh" #include "element_class_pentahedron_15_inline_impl.hh" #include "element_class_point_1_inline_impl.hh" #include "element_class_quadrangle_4_inline_impl.hh" #include "element_class_quadrangle_8_inline_impl.hh" #include "element_class_segment_2_inline_impl.hh" #include "element_class_segment_3_inline_impl.hh" #include "element_class_tetrahedron_10_inline_impl.hh" #include "element_class_tetrahedron_4_inline_impl.hh" #include "element_class_triangle_3_inline_impl.hh" #include "element_class_triangle_6_inline_impl.hh" /* -------------------------------------------------------------------------- */ #if defined(AKANTU_STRUCTURAL_MECHANICS) #include "element_class_structural.hh" #endif #if defined(AKANTU_COHESIVE_ELEMENT) #include "cohesive_element.hh" #endif #if defined(AKANTU_IGFEM) #include "element_class_igfem.hh" #endif #endif /* AKANTU_ELEMENT_CLASS_HH_ */ diff --git a/src/fe_engine/element_class_helper.hh b/src/fe_engine/element_class_helper.hh new file mode 100644 index 000000000..379a82433 --- /dev/null +++ b/src/fe_engine/element_class_helper.hh @@ -0,0 +1,75 @@ +#ifndef ELEMENT_CLASS_HELPER_HH +#define ELEMENT_CLASS_HELPER_HH + +#include "element_class.hh" + +namespace akantu { + +/* -------------------------------------------------------------------------- */ +template class ElementClassHelper {}; + +/* -------------------------------------------------------------------------- */ +template <> class ElementClassHelper<_ek_regular> { +public: + static inline Vector getN(const Vector & natural_coords, + ElementType type) { +#define GET_SHAPE_NATURAL(type) \ + auto nb_nodes_per_element = ElementClass::getNbNodesPerElement(); \ + Vector shapes(nb_nodes_per_element); \ + ElementClass::computeShapes(natural_coords, shapes); \ + return shapes + AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(GET_SHAPE_NATURAL); +#undef GET_SHAPE_NATURAL + + return Vector(0); + } + + /* ------------------------------------------------------------------------ */ + static inline Matrix getDNDS(const Vector & natural_coords, + ElementType type) { +#define GET_DNDS_NATURAL(type) \ + auto nb_nodes_per_element = ElementClass::getNbNodesPerElement(); \ + Matrix dnds(natural_coords.size(), nb_nodes_per_element); \ + ElementClass::computeDNDS(natural_coords, dnds); \ + return dnds + AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(GET_DNDS_NATURAL); +#undef GET_DNDS_NATURAL + + return Matrix(0, 0); + } + + /* ------------------------------------------------------------------------ */ + static inline Matrix getD2NDS2(const Vector & natural_coords, + ElementType type) { +#define GET_D2ND2S_NATURAL(type) \ + auto nb_nodes_per_element = ElementClass::getNbNodesPerElement(); \ + auto dim = natural_coords.size(); \ + Matrix d2nds2(dim * dim, nb_nodes_per_element); \ + ElementClass::computeD2NDS2(natural_coords, d2nds2); \ + return d2nds2 + AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(GET_D2ND2S_NATURAL); +#undef GET_D2NDS2_NATURAL + + return Matrix(0, 0); + } + + /* ------------------------------------------------------------------------ */ + static inline Matrix getJMat(const Vector & natural_coords, + const Matrix & positions, + ElementType type) { +#define GET_JMAT_NATURAL(type) \ + auto nb_nodes_per_element = ElementClass::getNbNodesPerElement(); \ + Matrix dnds(natural_coords.size(), nb_nodes_per_element); \ + Matrix jmat(dnds.rows(), positions.rows()); \ + ElementClass::computeDNDS(natural_coords, dnds); \ + ElementClass::computeJMat(dnds, positions, jmat); \ + return jmat + AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(GET_JMAT_NATURAL); +#undef GET_JMAT_NATURAL + + return Matrix(0, 0); + } +}; + +} // namespace akantu +#endif // ELEMENT_CLASS_HELPER_H diff --git a/src/fe_engine/element_class_tmpl.hh b/src/fe_engine/element_class_tmpl.hh index 89051d039..8fa0a9797 100644 --- a/src/fe_engine/element_class_tmpl.hh +++ b/src/fe_engine/element_class_tmpl.hh @@ -1,546 +1,542 @@ /** * @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: Wed Nov 29 2017 * * @brief Implementation of the inline templated function of the element class * descriptions * * * Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "element_class.hh" #include "gauss_integration_tmpl.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_ELEMENT_CLASS_TMPL_HH_ #define AKANTU_ELEMENT_CLASS_TMPL_HH_ namespace akantu { template inline constexpr auto ElementClass::getFacetTypes() { return VectorProxy( element_class_extra_geom_property::facet_type.data(), geometrical_element::getNbFacetTypes()); } /* -------------------------------------------------------------------------- */ /* GeometricalElement */ /* -------------------------------------------------------------------------- */ template inline constexpr auto GeometricalElement::getFacetLocalConnectivityPerElement(UInt t) { int pos = 0; for (UInt i = 0; i < t; ++i) { pos += geometrical_property::nb_facets[i] * geometrical_property::nb_nodes_per_facet[i]; } return MatrixProxy( geometrical_property::facet_connectivity_vect.data() + pos, geometrical_property::nb_facets[t], geometrical_property::nb_nodes_per_facet[t]); } /* -------------------------------------------------------------------------- */ template inline UInt GeometricalElement::getNbFacetsPerElement() { UInt total_nb_facets = 0; for (UInt n = 0; n < geometrical_property::nb_facet_types; ++n) { total_nb_facets += geometrical_property::nb_facets[n]; } return total_nb_facets; } /* -------------------------------------------------------------------------- */ template inline UInt GeometricalElement::getNbFacetsPerElement(UInt t) { return geometrical_property::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(), InterpolationProperty::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( InterpolationProperty::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( InterpolationProperty::natural_space_dimension, InterpolationProperty::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, UInt max_iterations, Real tolerance) { + Vector & natural_coords, UInt max_iterations, 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.zero(); // real space coordinates provided by initial guess Matrix physical_guess(spatial_dimension /*changed from - dimension */, 1); + dimension */ + , + 1); // objective function f = real_coords - physical_guess Matrix f(spatial_dimension /*changed from 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(dimension, spatial_dimension); // J^t Matrix Jt(spatial_dimension, dimension); - // G = J^t * J Matrix G(dimension, dimension); // Ginv = G^{-1} Matrix Ginv(dimension, dimension); // J = Ginv * J^t Matrix F(spatial_dimension, dimension); // dxi = \xi_{k+1} - \xi in the iterative process Matrix dxi(dimension, 1); Matrix dxit(1, dimension); - /* --------------------------- */ /* init before iteration loop */ /* --------------------------- */ // do interpolation auto update_f = [&f, &physical_guess, &natural_coords, &node_coords, - &mreal_coords, dimension, spatial_dimension]() { - Vector physical_guess_v(physical_guess.storage(), spatial_dimension /* changes - from dimension */); + &mreal_coords, spatial_dimension]() { + Vector physical_guess_v(physical_guess.storage(), spatial_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 auto error = f.norm(); return error; }; auto inverse_map_error = update_f(); - /* --------------------------- */ /* iteration loop */ /* --------------------------- */ UInt iterations{0}; while (tolerance < inverse_map_error and iterations < max_iterations) { // compute J^t interpolation_element::gradientOnNaturalCoordinates(natural_coords, node_coords, Jt); J = Jt.transpose(); - + // compute G G.mul(J, J); // inverse G Ginv.inverse(G); // compute F F.mul(J, Ginv); // compute increment dxit.mul(f, F); dxi = dxit.transpose(); - + // update our guess natural_coords += Vector(dxi(0)); inverse_map_error = update_f(); iterations++; } - // 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, UInt max_iterations, 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, max_iterations, tolerance); } } } // namespace akantu #endif /* AKANTU_ELEMENT_CLASS_TMPL_HH_ */ diff --git a/src/fe_engine/element_classes/element_class_quadrangle_4_inline_impl.hh b/src/fe_engine/element_classes/element_class_quadrangle_4_inline_impl.hh index 015be3bd6..f2468d706 100644 --- a/src/fe_engine/element_classes/element_class_quadrangle_4_inline_impl.hh +++ b/src/fe_engine/element_classes/element_class_quadrangle_4_inline_impl.hh @@ -1,182 +1,182 @@ /** * @file element_class_quadrangle_4_inline_impl.hh * * @author Nicolas Richart * * @date creation: Mon Dec 13 2010 * @date last modification: Wed Oct 11 2017 * * @brief Specialization of the element_class class for the type _quadrangle_4 * * * Copyright (©) 2010-2018 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 . * * * @verbatim \eta ^ (-1,1) | (1,1) x---------x | | | | | | --|---------|-----> \xi | | | | | | x---------x (-1,-1) | (1,-1) @endverbatim * * @f[ * \begin{array}{lll} * N1 = (1 - \xi) (1 - \eta) / 4 * & \frac{\partial N1}{\partial \xi} = - (1 - \eta) / 4 * & \frac{\partial N1}{\partial \eta} = - (1 - \xi) / 4 \\ * N2 = (1 + \xi) (1 - \eta) / 4 \\ * & \frac{\partial N2}{\partial \xi} = (1 - \eta) / 4 * & \frac{\partial N2}{\partial \eta} = - (1 + \xi) / 4 \\ * N3 = (1 + \xi) (1 + \eta) / 4 \\ * & \frac{\partial N3}{\partial \xi} = (1 + \eta) / 4 * & \frac{\partial N3}{\partial \eta} = (1 + \xi) / 4 \\ * N4 = (1 - \xi) (1 + \eta) / 4 * & \frac{\partial N4}{\partial \xi} = - (1 + \eta) / 4 * & \frac{\partial N4}{\partial \eta} = (1 - \xi) / 4 \\ * \end{array} * @f] * * @f{eqnarray*}{ * \xi_{q0} &=& 0 \qquad \eta_{q0} = 0 * @f} */ /* -------------------------------------------------------------------------- */ #include "element_class.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ AKANTU_DEFINE_ELEMENT_CLASS_PROPERTY(_quadrangle_4, _gt_quadrangle_4, _itp_lagrange_quadrangle_4, _ek_regular, 2, _git_segment, 2); /* -------------------------------------------------------------------------- */ template <> template inline void InterpolationElement<_itp_lagrange_quadrangle_4>::computeShapes( const vector_type & c, vector_type & N) { N(0) = 1. / 4. * (1. - c(0)) * (1. - c(1)); /// N1(q_0) N(1) = 1. / 4. * (1. + c(0)) * (1. - c(1)); /// N2(q_0) N(2) = 1. / 4. * (1. + c(0)) * (1. + c(1)); /// N3(q_0) N(3) = 1. / 4. * (1. - c(0)) * (1. + c(1)); /// N4(q_0) } /* -------------------------------------------------------------------------- */ template <> template inline void InterpolationElement<_itp_lagrange_quadrangle_4>::computeDNDS( const vector_type & c, matrix_type & dnds) { /** * @f[ * dnds = \left( * \begin{array}{cccc} * \frac{\partial N1}{\partial \xi} & \frac{\partial N2}{\partial * \xi} * & \frac{\partial N3}{\partial \xi} & \frac{\partial * N4}{\partial \xi}\\ * \frac{\partial N1}{\partial \eta} & \frac{\partial N2}{\partial * \eta} * & \frac{\partial N3}{\partial \eta} & \frac{\partial * N4}{\partial \eta} * \end{array} * \right) * @f] */ dnds(0, 0) = -1. / 4. * (1. - c(1)); dnds(0, 1) = 1. / 4. * (1. - c(1)); dnds(0, 2) = 1. / 4. * (1. + c(1)); dnds(0, 3) = -1. / 4. * (1. + c(1)); dnds(1, 0) = -1. / 4. * (1. - c(0)); dnds(1, 1) = -1. / 4. * (1. + c(0)); dnds(1, 2) = 1. / 4. * (1. + c(0)); dnds(1, 3) = 1. / 4. * (1. - c(0)); } /* -------------------------------------------------------------------------- */ template<> template -inline void InterpolationElement<_itp_lagrange_quadrangle_4>::computeDN2DS2( - const vector_type & /*c*/, matrix_type & dn2ds2) { +inline void InterpolationElement<_itp_lagrange_quadrangle_4>::computeD2NDS2( + const vector_type & /*c*/, matrix_type & d2nds2) { + d2nds2.zero(); + // d2nds2(0, 0) = 0; + // d2nds2(0, 1) = 0; + // d2nds2(0, 2) = 0; + // d2nds2(0, 3) = 0; - dn2ds2(0, 0) = 0; - dn2ds2(0, 1) = 0; - dn2ds2(0, 2) = 0; - dn2ds2(0, 3) = 0; + d2nds2(1, 0) = 1./4.; + d2nds2(1, 1) = -1./4.; + d2nds2(1, 2) = 1./4.; + d2nds2(1, 3) = -1./4.; - dn2ds2(1, 0) = 1./4.; - dn2ds2(1, 1) = -1./4.; - dn2ds2(1, 2) = 1./4.; - dn2ds2(1, 3) = -1./4.; + d2nds2(2, 0) = 1./4.; + d2nds2(2, 1) = -1./4.; + d2nds2(2, 2) = 1./4.; + d2nds2(2, 3) = -1./4.; - dn2ds2(2, 0) = 1./4.; - dn2ds2(2, 1) = -1./4.; - dn2ds2(2, 2) = 1./4.; - dn2ds2(2, 3) = -1./4.; - - dn2ds2(3, 0) = 0; - dn2ds2(3, 1) = 0; - dn2ds2(3, 2) = 0; - dn2ds2(3, 3) = 0; + // d2nds2(3, 0) = 0; + // d2nds2(3, 1) = 0; + // d2nds2(3, 2) = 0; + // d2nds2(3, 3) = 0; } /* -------------------------------------------------------------------------- */ template <> inline void InterpolationElement<_itp_lagrange_quadrangle_4>::computeSpecialJacobian( const Matrix & J, Real & jac) { Vector vprod(J.cols()); Matrix Jt(J.transpose(), true); vprod.crossProduct(Jt(0), Jt(1)); jac = vprod.norm(); } /* -------------------------------------------------------------------------- */ template <> inline Real GeometricalElement<_gt_quadrangle_4>::getInradius(const Matrix & coord) { Vector u0 = coord(0); Vector u1 = coord(1); Vector u2 = coord(2); Vector u3 = coord(3); Real a = u0.distance(u1); Real b = u1.distance(u2); Real c = u2.distance(u3); Real d = u3.distance(u0); // Real septimetre = (a + b + c + d) / 2.; // Real p = Math::distance_2d(coord + 0, coord + 4); // Real q = Math::distance_2d(coord + 2, coord + 6); // Real area = sqrt(4*(p*p * q*q) - (a*a + b*b + c*c + d*d)*(a*a + c*c - b*b - // d*d)) / 4.; // Real h = sqrt(area); // to get a length // Real h = area / septimetre; // formula of inradius for circumscritable // quadrelateral Real h = std::min({a, b, c, d}); return h; } } // namespace akantu diff --git a/src/fe_engine/element_classes/element_class_segment_2_inline_impl.hh b/src/fe_engine/element_classes/element_class_segment_2_inline_impl.hh index d6e61bc21..0f9293a78 100644 --- a/src/fe_engine/element_classes/element_class_segment_2_inline_impl.hh +++ b/src/fe_engine/element_classes/element_class_segment_2_inline_impl.hh @@ -1,128 +1,120 @@ /** * @file element_class_segment_2_inline_impl.hh * * @author Nicolas Richart * * @date creation: Fri Jul 16 2010 * @date last modification: Wed Oct 11 2017 * * @brief Specialization of the element_class class for the type _segment_2 * * * Copyright (©) 2010-2018 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 . * * * @verbatim q --x--------|--------x---> x -1 0 1 @endverbatim * * @f{eqnarray*}{ * w_1(x) &=& 1/2(1 - x) \\ * w_2(x) &=& 1/2(1 + x) * @f} * * @f{eqnarray*}{ * x_{q} &=& 0 * @f} */ /* -------------------------------------------------------------------------- */ #include "element_class.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ AKANTU_DEFINE_ELEMENT_CLASS_PROPERTY(_segment_2, _gt_segment_2, _itp_lagrange_segment_2, _ek_regular, 1, _git_segment, 1); /* -------------------------------------------------------------------------- */ template <> template inline void InterpolationElement<_itp_lagrange_segment_2>::computeShapes( const vector_type & natural_coords, vector_type & N) { /// natural coordinate Real c = natural_coords(0); /// shape functions N(0) = 0.5 * (1 - c); N(1) = 0.5 * (1 + c); } /* -------------------------------------------------------------------------- */ template <> template inline void InterpolationElement<_itp_lagrange_segment_2>::computeDNDS( __attribute__((unused)) const vector_type & natural_coords, matrix_type & dnds) { /// dN1/de dnds(0, 0) = -.5; /// dN2/de dnds(0, 1) = .5; } /* -------------------------------------------------------------------------- */ template<> template -inline void InterpolationElement<_itp_lagrange_segment_2>::computeDN2DS2( - __attribute__((unused)) const vector_type & natural_coords, matrix_type & dn2ds2) { - - /** - * @f[ - * dn2ds2 = \left( - * \begin{array}{} - * @f] - */ - - dn2ds2 *= 0; +inline void InterpolationElement<_itp_lagrange_segment_2>::computeD2NDS2( + const vector_type & /*natural_coords*/, matrix_type & d2nds2) { + d2nds2.zero(); } /* -------------------------------------------------------------------------- */ template <> inline void InterpolationElement<_itp_lagrange_segment_2>::computeSpecialJacobian( const Matrix & dxds, Real & jac) { jac = dxds.norm(); } /* -------------------------------------------------------------------------- */ template <> inline Real GeometricalElement<_gt_segment_2>::getInradius(const Matrix & coord) { if (coord.rows() == 2) { return Math::distance_2d(coord(0).storage(), coord(1).storage()); } if (coord.rows() == 3) { return Math::distance_3d(coord(0).storage(), coord(1).storage()); } return std::abs(coord(0, 0) - coord(0, 1)); } // /* -------------------------------------------------------------------------- // */ // template<> inline bool ElementClass<_segment_2>::contains(const Vector // & natural_coords) { // if (natural_coords(0) < -1.) return false; // if (natural_coords(0) > 1.) return false; // return true; // } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/fe_engine/element_classes/element_class_triangle_3_inline_impl.hh b/src/fe_engine/element_classes/element_class_triangle_3_inline_impl.hh index 9f7546ed6..2a632eb54 100644 --- a/src/fe_engine/element_classes/element_class_triangle_3_inline_impl.hh +++ b/src/fe_engine/element_classes/element_class_triangle_3_inline_impl.hh @@ -1,154 +1,146 @@ /** * @file element_class_triangle_3_inline_impl.hh * * @author Nicolas Richart * * @date creation: Fri Jul 16 2010 * @date last modification: Wed Oct 11 2017 * * @brief Specialization of the element_class class for the type _triangle_3 * * * Copyright (©) 2010-2018 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 . * * * @verbatim \eta ^ | x (0,0,1) |` | ` | q ` | ° ` x--------x-----> \xi (1,0,0) (0,1,0) @endverbatim * * @f{eqnarray*}{ * N1 &=& 1 - \xi - \eta \\ * N2 &=& \xi \\ * N3 &=& \eta * @f} * * @f{eqnarray*}{ * \xi_{q0} &=& 1/3 \qquad \eta_{q0} = 1/3 * @f} */ /* -------------------------------------------------------------------------- */ #include "element_class.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ AKANTU_DEFINE_ELEMENT_CLASS_PROPERTY(_triangle_3, _gt_triangle_3, _itp_lagrange_triangle_3, _ek_regular, 2, _git_triangle, 1); /* -------------------------------------------------------------------------- */ template <> template inline void InterpolationElement<_itp_lagrange_triangle_3>::computeShapes( const vector_type & natural_coords, vector_type & N) { /// Natural coordinates Real c0 = 1 - natural_coords(0) - natural_coords(1); /// @f$ c0 = 1 - \xi - \eta @f$ Real c1 = natural_coords(0); /// @f$ c1 = \xi @f$ Real c2 = natural_coords(1); /// @f$ c2 = \eta @f$ N(0) = c0; /// N1(q_0) N(1) = c1; /// N2(q_0) N(2) = c2; /// N3(q_0) } /* -------------------------------------------------------------------------- */ template <> template inline void InterpolationElement<_itp_lagrange_triangle_3>::computeDNDS( __attribute__((unused)) const vector_type & natural_coords, matrix_type & dnds) { /** * @f[ * dnds = \left( * \begin{array}{cccccc} * \frac{\partial N1}{\partial \xi} & \frac{\partial N2}{\partial * \xi} & \frac{\partial N3}{\partial \xi} \\ * \frac{\partial N1}{\partial \eta} & \frac{\partial N2}{\partial * \eta} & \frac{\partial N3}{\partial \eta} * \end{array} * \right) * @f] */ dnds(0, 0) = -1.; dnds(0, 1) = 1.; dnds(0, 2) = 0.; dnds(1, 0) = -1.; dnds(1, 1) = 0.; dnds(1, 2) = 1.; } /* -------------------------------------------------------------------------- */ -template<> +template <> template -inline void InterpolationElement<_itp_lagrange_triangle_3>::computeDN2DS2( - __attribute__((unused)) const vector_type & natural_coords, matrix_type & dn2ds2) { - - /** - * @f[ - * dn2ds2 = \left( - * \begin{array}{} - * @f] - */ - - dn2ds2 *= 0; +inline void InterpolationElement<_itp_lagrange_triangle_3>::computeD2NDS2( + const vector_type & /*natural_coords*/, matrix_type & d2nds2) { + d2nds2.zero(); } /* -------------------------------------------------------------------------- */ template <> inline void InterpolationElement<_itp_lagrange_triangle_3>::computeSpecialJacobian( const Matrix & J, Real & jac) { Vector vprod(J.cols()); Matrix Jt(J.transpose(), true); vprod.crossProduct(Jt(0), Jt(1)); jac = vprod.norm(); } /* -------------------------------------------------------------------------- */ template <> inline Real GeometricalElement<_gt_triangle_3>::getInradius(const Matrix & coord) { if (coord.rows() == 2) { - return 2. * Math::triangle_inradius_2d(coord(0).storage(), coord(1).storage(), - coord(2).storage()); + return 2. * Math::triangle_inradius_2d( + coord(0).storage(), coord(1).storage(), coord(2).storage()); } return 2. * Math::triangle_inradius_3d(coord(0).storage(), coord(1).storage(), - coord(2).storage()); + coord(2).storage()); } /* -------------------------------------------------------------------------- */ // template<> inline bool ElementClass<_triangle_3>::contains(const Vector // & natural_coords) { // if (natural_coords[0] < 0.) return false; // if (natural_coords[0] > 1.) return false; // if (natural_coords[1] < 0.) return false; // if (natural_coords[1] > 1.) return false; // if (natural_coords[0]+natural_coords[1] > 1.) return false; // return true; // } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/fe_engine/fe_engine_inline_impl.hh b/src/fe_engine/fe_engine_inline_impl.hh index 4e0e57288..87968ecaf 100644 --- a/src/fe_engine/fe_engine_inline_impl.hh +++ b/src/fe_engine/fe_engine_inline_impl.hh @@ -1,194 +1,190 @@ /** * @file fe_engine_inline_impl.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Tue Jul 20 2010 * @date last modification: Sun Aug 13 2017 * * @brief Implementation of the inline functions of the FEEngine Class * * * Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "element_class.hh" #include "fe_engine.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ #include "element_type_conversion.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_FE_ENGINE_INLINE_IMPL_HH_ #define AKANTU_FE_ENGINE_INLINE_IMPL_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ inline Real FEEngine::getElementInradius(const Matrix & coord, ElementType type) { Real inradius = 0; #define GET_INRADIUS(type) inradius = ElementClass::getInradius(coord); AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_INRADIUS); #undef GET_INRADIUS return inradius; } /* -------------------------------------------------------------------------- */ -inline InterpolationType -FEEngine::getInterpolationType(ElementType type) { +inline InterpolationType FEEngine::getInterpolationType(ElementType type) { return convertType(type); } /* -------------------------------------------------------------------------- */ /// @todo rewrite this function in order to get the cohesive element /// type directly from the facet #if defined(AKANTU_COHESIVE_ELEMENT) inline ElementType FEEngine::getCohesiveElementType(ElementType type) { ElementType ctype; #define GET_COHESIVE_TYPE(type) \ ctype = CohesiveFacetProperty::cohesive_type; AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_COHESIVE_TYPE); #undef GET_COHESIVE_TYPE return ctype; } #else -inline ElementType -FEEngine::getCohesiveElementType(__attribute__((unused)) - ElementType type_facet) { +inline ElementType FEEngine::getCohesiveElementType(__attribute__((unused)) + ElementType type_facet) { return _not_defined; } #endif /* -------------------------------------------------------------------------- */ #if defined(AKANTU_IGFEM) } // akantu #include "igfem_helper.hh" namespace akantu { -inline Vector -FEEngine::getIGFEMElementTypes(ElementType type) { +inline Vector FEEngine::getIGFEMElementTypes(ElementType type) { #define GET_IGFEM_ELEMENT_TYPES(type) \ return IGFEMHelper::getIGFEMElementTypes(); AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(GET_IGFEM_ELEMENT_TYPES); #undef GET_IGFEM_ELEMENT_TYPES } #endif /* -------------------------------------------------------------------------- */ template void FEEngine::extractNodalToElementField(const Mesh & mesh, const Array & nodal_f, Array & elemental_f, ElementType type, GhostType ghost_type, const Array & filter_elements) { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_degree_of_freedom = nodal_f.getNbComponent(); UInt nb_element = mesh.getNbElement(type, ghost_type); UInt * conn_val = mesh.getConnectivity(type, ghost_type).storage(); if (filter_elements != empty_filter) { nb_element = filter_elements.size(); } elemental_f.resize(nb_element); T * nodal_f_val = nodal_f.storage(); T * f_val = elemental_f.storage(); UInt * el_conn; for (UInt el = 0; el < nb_element; ++el) { if (filter_elements != empty_filter) { el_conn = conn_val + filter_elements(el) * nb_nodes_per_element; } else { el_conn = conn_val + el * nb_nodes_per_element; } for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = *(el_conn + n); std::copy(nodal_f_val + node * nb_degree_of_freedom, nodal_f_val + (node + 1) * nb_degree_of_freedom, f_val); f_val += nb_degree_of_freedom; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void FEEngine::filterElementalData(const Mesh & mesh, const Array & elem_f, - Array & filtered_f, - ElementType type, + Array & filtered_f, ElementType type, GhostType ghost_type, const Array & filter_elements) { AKANTU_DEBUG_IN(); UInt nb_element = mesh.getNbElement(type, ghost_type); if (nb_element == 0) { filtered_f.resize(0); return; } UInt nb_degree_of_freedom = elem_f.getNbComponent(); UInt nb_data_per_element = elem_f.size() / nb_element; if (filter_elements != empty_filter) { nb_element = filter_elements.size(); } filtered_f.resize(nb_element * nb_data_per_element); T * elem_f_val = elem_f.storage(); T * f_val = filtered_f.storage(); UInt el_offset; for (UInt el = 0; el < nb_element; ++el) { if (filter_elements != empty_filter) { el_offset = filter_elements(el); } else { el_offset = el; } std::copy(elem_f_val + el_offset * nb_data_per_element * nb_degree_of_freedom, elem_f_val + (el_offset + 1) * nb_data_per_element * nb_degree_of_freedom, f_val); f_val += nb_degree_of_freedom * nb_data_per_element; } AKANTU_DEBUG_OUT(); } } // namespace akantu #endif /* AKANTU_FE_ENGINE_INLINE_IMPL_HH_ */ diff --git a/src/model/contact_mechanics/contact_detector.cc b/src/model/contact_mechanics/contact_detector.cc index 25a89c55d..20b8eb9cb 100644 --- a/src/model/contact_mechanics/contact_detector.cc +++ b/src/model/contact_mechanics/contact_detector.cc @@ -1,362 +1,353 @@ /** * @file contact_detector.cc * * @author Mohit Pundir * * @date creation: Wed Sep 12 2018 * @date last modification: Fri Sep 21 2018 * * @brief Mother class for all detection algorithms * * @section LICENSE * * Copyright (©) 2010-2018 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 "contact_detector.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ ContactDetector::ContactDetector(Mesh & mesh, const ID & id, const MemoryID & memory_id) : ContactDetector(mesh, mesh.getNodes(), id, memory_id) {} /* -------------------------------------------------------------------------- */ ContactDetector::ContactDetector(Mesh & mesh, Array positions, const ID & id, const MemoryID & memory_id) : Memory(id, memory_id), Parsable(ParserType::_contact_detector, id), mesh(mesh), positions(0, mesh.getSpatialDimension()) { AKANTU_DEBUG_IN(); this->spatial_dimension = mesh.getSpatialDimension(); this->positions.copy(positions); this->parseSection(); AKANTU_DEBUG_OUT(); } - /* -------------------------------------------------------------------------- */ void ContactDetector::parseSection() { const Parser & parser = getStaticParser(); const ParserSection & section = *(parser.getSubSections(ParserType::_contact_detector).first); auto type = section.getParameterValue("type"); if (type == "implicit") { this->detection_type = _implicit; } else if (type == "explicit") { this->detection_type = _explicit; - } else { + } else { AKANTU_ERROR("Unknown detection type : " << type); } this->projection_tolerance = - section.getParameterValue("projection_tolerance"); - this->max_iterations = - section.getParameterValue("max_iterations"); + section.getParameterValue("projection_tolerance"); + this->max_iterations = section.getParameterValue("max_iterations"); this->extension_tolerance = - section.getParameterValue("extension_tolerance"); + section.getParameterValue("extension_tolerance"); } /* -------------------------------------------------------------------------- */ void ContactDetector::search(Array & elements, - Array & gaps, Array & normals, - Array & tangents, - Array & projections) { + Array & gaps, Array & normals, + Array & tangents, + Array & projections) { UInt surface_dimension = spatial_dimension - 1; this->mesh.fillNodesToElements(surface_dimension); this->computeMaximalDetectionDistance(); contact_pairs.clear(); SpatialGrid master_grid(spatial_dimension); SpatialGrid slave_grid(spatial_dimension); this->globalSearch(slave_grid, master_grid); this->localSearch(slave_grid, master_grid); this->createContactElements(elements, gaps, normals, tangents, projections); - } - /* -------------------------------------------------------------------------- */ void ContactDetector::globalSearch(SpatialGrid & slave_grid, SpatialGrid & master_grid) { auto & master_list = surface_selector->getMasterList(); auto & slave_list = surface_selector->getSlaveList(); BBox bbox_master(spatial_dimension); this->constructBoundingBox(bbox_master, master_list); BBox bbox_slave(spatial_dimension); this->constructBoundingBox(bbox_slave, slave_list); auto && bbox_intersection = bbox_master.intersection(bbox_slave); AKANTU_DEBUG_INFO("Intersection BBox " << bbox_intersection); Vector center(spatial_dimension); bbox_intersection.getCenter(center); Vector spacing(spatial_dimension); this->computeCellSpacing(spacing); master_grid.setCenter(center); master_grid.setSpacing(spacing); this->constructGrid(master_grid, bbox_intersection, master_list); slave_grid.setCenter(center); slave_grid.setSpacing(spacing); this->constructGrid(slave_grid, bbox_intersection, slave_list); // search slave grid nodes in contactelement array and if they exits // and still have orthogonal projection on its associated master // facetremove it from the spatial grid or do not consider it for // local search, maybe better option will be to have spatial grid of // type node info and one of the variable of node info should be // facet already exits // so contact elements will be updated based on the above // consideration , this means only those contact elements will be // keep whose slave node is still in intersection bbox and still has // projection in its master facet // also if slave node is already exists in contact element and // orthogonal projection does not exits then search the associated // master facets with the current master facets within a given // radius , this is subjected to computational cost as searching // neighbbor cells can be more effective. } /* -------------------------------------------------------------------------- */ void ContactDetector::localSearch(SpatialGrid & slave_grid, SpatialGrid & master_grid) { // local search // out of these array check each cell for closet node in that cell // and neighbouring cells find the actual orthogonally closet // check the projection of slave node on master facets connected to // the closet master node, if yes update the contact element with // slave node and master node and master surfaces connected to the // master node // these master surfaces will be needed later to update contact // elements /// find the closet master node for each slave node for (auto && cell_id : slave_grid) { /// loop over all the slave nodes of the current cell for (auto && slave_node : slave_grid.getCell(cell_id)) { bool pair_exists = false; Vector pos(spatial_dimension); for (UInt s : arange(spatial_dimension)) pos(s) = this->positions(slave_node, s); Real closet_distance = std::numeric_limits::max(); UInt closet_master_node; /// loop over all the neighboring cells of the current cell for (auto && neighbor_cell : cell_id.neighbors()) { /// loop over the data of neighboring cells from master grid for (auto && master_node : master_grid.getCell(neighbor_cell)) { /// check for self contact if (slave_node == master_node) continue; - bool is_valid = true; + bool is_valid = true; Array elements; this->mesh.getAssociatedElements(slave_node, elements); for (auto & elem : elements) { if (elem.kind() != _ek_regular) continue; Vector connectivity = const_cast(this->mesh).getConnectivity(elem); auto node_iter = std::find(connectivity.begin(), connectivity.end(), master_node); if (node_iter != connectivity.end()) { is_valid = false; - break; + break; } } Vector pos2(spatial_dimension); for (UInt s : arange(spatial_dimension)) pos2(s) = this->positions(master_node, s); Real distance = pos.distance(pos2); if (distance <= closet_distance and is_valid) { closet_master_node = master_node; closet_distance = distance; pair_exists = true; } } } if (pair_exists) contact_pairs.push_back(std::make_pair(slave_node, closet_master_node)); } } } - - /* -------------------------------------------------------------------------- */ -void ContactDetector::createContactElements(Array & contact_elements, - Array & gaps, Array & normals, - Array & tangents, Array & projections) { +void ContactDetector::createContactElements( + Array & contact_elements, Array & gaps, + Array & normals, Array & tangents, Array & projections) { auto surface_dimension = spatial_dimension - 1; - + Real alpha; switch (detection_type) { case _explicit: { alpha = 1.0; break; } case _implicit: { alpha = -1.0; break; - } + } default: - AKANTU_EXCEPTION(detection_type << " is not a valid contact detection type"); + AKANTU_EXCEPTION(detection_type + << " is not a valid contact detection type"); break; } - + for (auto & pairs : contact_pairs) { const auto & slave_node = pairs.first; Vector slave(spatial_dimension); for (UInt s : arange(spatial_dimension)) slave(s) = this->positions(slave_node, s); const auto & master_node = pairs.second; Array elements; this->mesh.getAssociatedElements(master_node, elements); auto & gap = gaps.begin()[slave_node]; Vector normal(normals.begin(spatial_dimension)[slave_node]); Vector projection(projections.begin(surface_dimension)[slave_node]); - Matrix tangent(tangents.begin(surface_dimension, - spatial_dimension)[slave_node]); - auto index = GeometryUtils::orthogonalProjection(mesh, positions, slave, elements, - gap, projection, normal, - tangent, - alpha, - this->max_iterations, - this->projection_tolerance, - this->extension_tolerance); + Matrix tangent( + tangents.begin(surface_dimension, spatial_dimension)[slave_node]); + auto index = GeometryUtils::orthogonalProjection( + mesh, positions, slave, elements, gap, projection, normal, tangent, + alpha, this->max_iterations, this->projection_tolerance, + this->extension_tolerance); // if a valid projection is not found on the patch of elements // index is -1 or if not a valid self contact, the contact element // is not created - if (index == UInt(-1) or !isValidSelfContact(slave_node, gap, normal)) { + if (index == UInt(-1) or !isValidSelfContact(slave_node, gap, normal)) { gap *= 0; normal *= 0; projection *= 0; tangent *= 0; continue; } - + // create contact element contact_elements.push_back(ContactElement(slave_node, elements[index])); } contact_pairs.clear(); } /* -------------------------------------------------------------------------- */ -/*void ContactDetector::createContactElements(Array & contact_elements, - Array & gaps, Array & normals, - Array & projections) { +/*void ContactDetector::createContactElements(Array & + contact_elements, Array & gaps, Array & normals, Array & + projections) { auto surface_dimension = spatial_dimension - 1; - + Real alpha; switch (detection_type) { case _explicit: { alpha = 1.0; break; } case _implicit: { alpha = -1.0; break; - } + } default: - AKANTU_EXCEPTION(detection_type << " is not a valid contact detection type"); - break; + AKANTU_EXCEPTION(detection_type << " is not a valid contact detection + type"); break; } - + for (auto & pairs : contact_pairs) { const auto & slave_node = pairs.first; Vector slave(spatial_dimension); for (UInt s : arange(spatial_dimension)) slave(s) = this->positions(slave_node, s); const auto & master_node = pairs.second; Array elements; this->mesh.getAssociatedElements(master_node, elements); - + auto & gap = gaps.begin()[slave_node]; Vector normal(normals.begin(spatial_dimension)[slave_node]); Vector projection(projections.begin(surface_dimension)[slave_node]); - auto index = GeometryUtils::orthogonalProjection(mesh, positions, slave, elements, - gap, projection, normal, alpha, - this->max_iterations, - this->projection_tolerance, - this->extension_tolerance); + auto index = GeometryUtils::orthogonalProjection(mesh, positions, slave, + elements, gap, projection, normal, alpha, this->max_iterations, + this->projection_tolerance, + this->extension_tolerance); // if a valid projection is not found on the patch of elements // index is -1 or if not a valid self contact, the contact element // is not created - if (index == UInt(-1) or !isValidSelfContact(slave_node, gap, normal)) { + if (index == UInt(-1) or !isValidSelfContact(slave_node, gap, normal)) { gap *= 0; normal *= 0; continue; } - + // create contact element contact_elements.push_back(ContactElement(slave_node, elements[index])); } contact_pairs.clear(); }*/ } // namespace akantu diff --git a/src/model/contact_mechanics/geometry_utils.cc b/src/model/contact_mechanics/geometry_utils.cc index 7a66183bc..a41ec9496 100644 --- a/src/model/contact_mechanics/geometry_utils.cc +++ b/src/model/contact_mechanics/geometry_utils.cc @@ -1,790 +1,750 @@ /** * @file geometry_utils.cc * * @author Mohit Pundir * * @date creation: Mon Mmay 20 2019 * @date last modification: Mon May 20 2019 * * @brief Implementation of various utilities needed for contact geometry * * @section LICENSE * * Copyright (©) 2010-2018 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 "geometry_utils.hh" +#include "element_class_helper.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ void GeometryUtils::normal(const Mesh & mesh, const Array & positions, const Element & element, Vector & normal, bool outward) { UInt spatial_dimension = mesh.getSpatialDimension(); UInt surface_dimension = spatial_dimension - 1; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); Matrix coords(spatial_dimension, nb_nodes_per_element); UInt * elem_val = mesh.getConnectivity(element.type, _not_ghost).storage(); mesh.extractNodalValuesFromElement(positions, coords.storage(), elem_val + element.element * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); Matrix vectors(spatial_dimension, surface_dimension); switch (spatial_dimension) { case 1: { normal[0] = 1; break; } case 2: { vectors(0) = Vector(coords(1)) - Vector(coords(0)); Math::normal2(vectors.storage(), normal.storage()); break; } case 3: { vectors(0) = Vector(coords(1)) - Vector(coords(0)); vectors(1) = Vector(coords(2)) - Vector(coords(0)); Math::normal3(vectors(0).storage(), vectors(1).storage(), normal.storage()); break; } default: { AKANTU_ERROR("Unknown dimension : " << spatial_dimension); } } // to ensure that normal is always outwards from master surface if (outward) { const auto & element_to_subelement = mesh.getElementToSubelement(element.type)(element.element); Vector outside(spatial_dimension); mesh.getBarycenter(element, outside); // check if mesh facets exists for cohesive elements contact Vector inside(spatial_dimension); if (mesh.isMeshFacets()) { mesh.getMeshParent().getBarycenter(element_to_subelement[0], inside); } else { mesh.getBarycenter(element_to_subelement[0], inside); } Vector inside_to_outside = outside - inside; auto projection = inside_to_outside.dot(normal); if (projection < 0) { normal *= -1.0; } } } /* -------------------------------------------------------------------------- */ void GeometryUtils::normal(const Mesh & mesh, const Element & element, Matrix & tangents, Vector & normal, bool outward) { UInt spatial_dimension = mesh.getSpatialDimension(); // to ensure that normal is always outwards from master surface we // compute a direction vector form inside of element attached to the // suurface element Vector inside_to_outside(spatial_dimension); if (outward) { const auto & element_to_subelement = mesh.getElementToSubelement(element.type)(element.element); Vector outside(spatial_dimension); mesh.getBarycenter(element, outside); // check if mesh facets exists for cohesive elements contact Vector inside(spatial_dimension); if (mesh.isMeshFacets()) { mesh.getMeshParent().getBarycenter(element_to_subelement[0], inside); } else { mesh.getBarycenter(element_to_subelement[0], inside); } inside_to_outside = outside - inside; // auto projection = inside_to_outside.dot(normal); // if (projection < 0) { // normal *= -1.0; //} } // to ensure that direction of tangents are correct, cross product // of tangents should give be in the same direction as of inside to outside switch (spatial_dimension) { case 2: { normal[0] = -tangents(0, 1); normal[1] = tangents(0, 0); auto ddot = inside_to_outside.dot(normal); if (ddot < 0) { tangents *= -1.0; normal *= -1.0; } break; } case 3: { auto tang_trans = tangents.transpose(); auto tang1 = Vector(tang_trans(0)); auto tang2 = Vector(tang_trans(1)); auto tang1_cross_tang2 = tang1.crossProduct(tang2); normal = tang1_cross_tang2 / tang1_cross_tang2.norm(); auto ddot = inside_to_outside.dot(normal); if (ddot < 0) { tang_trans(1) *= -1.0; normal *= -1.0; } tangents = tang_trans.transpose(); break; } default: break; } } /* -------------------------------------------------------------------------- */ void GeometryUtils::covariantBasis(const Mesh & mesh, const Array & positions, const Element & element, const Vector & normal, Vector & natural_coord, Matrix & tangents) { UInt spatial_dimension = mesh.getSpatialDimension(); - const ElementType & type = element.type; + const ElementType type = element.type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); UInt * elem_val = mesh.getConnectivity(type, _not_ghost).storage(); Matrix nodes_coord(spatial_dimension, nb_nodes_per_element); mesh.extractNodalValuesFromElement(positions, nodes_coord.storage(), elem_val + element.element * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); - UInt surface_dimension = spatial_dimension - 1; - Matrix dnds(surface_dimension, nb_nodes_per_element); - -#define GET_SHAPE_DERIVATIVES_NATURAL(type) \ - ElementClass::computeDNDS(natural_coord, dnds) - AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); -#undef GET_SHAPE_DERIVATIVES_NATURAL - + auto && dnds = ElementClassHelper<_ek_regular>::getDNDS(natural_coord, type); tangents.mul(dnds, nodes_coord); auto temp_tangents = tangents.transpose(); for (UInt i = 0; i < spatial_dimension - 1; ++i) { auto temp = Vector(temp_tangents(i)); temp_tangents(i) = temp.normalize(); } tangents = temp_tangents.transpose(); // to ensure that direction of tangents are correct, cross product // of tangents should give the normal vector computed earlier switch (spatial_dimension) { case 2: { Vector e_z(3); e_z[0] = 0.; e_z[1] = 0.; e_z[2] = 1.; Vector tangent(3); tangent[0] = tangents(0, 0); tangent[1] = tangents(0, 1); tangent[2] = 0.; auto exp_normal = e_z.crossProduct(tangent); auto ddot = normal.dot(exp_normal); if (ddot < 0) { tangents *= -1.0; } break; } case 3: { auto tang_trans = tangents.transpose(); auto tang1 = Vector(tang_trans(0)); auto tang2 = Vector(tang_trans(1)); auto tang1_cross_tang2 = tang1.crossProduct(tang2); auto exp_normal = tang1_cross_tang2 / tang1_cross_tang2.norm(); auto ddot = normal.dot(exp_normal); if (ddot < 0) { tang_trans(1) *= -1.0; } tangents = tang_trans.transpose(); break; } default: break; } } /* -------------------------------------------------------------------------- */ void GeometryUtils::covariantBasis(const Mesh & mesh, const Array & positions, const Element & element, Vector & natural_coord, Matrix & tangents) { UInt spatial_dimension = mesh.getSpatialDimension(); const ElementType & type = element.type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); UInt * elem_val = mesh.getConnectivity(type, _not_ghost).storage(); Matrix nodes_coord(spatial_dimension, nb_nodes_per_element); mesh.extractNodalValuesFromElement(positions, nodes_coord.storage(), elem_val + element.element * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); - UInt surface_dimension = spatial_dimension - 1; - Matrix dnds(surface_dimension, nb_nodes_per_element); - -#define GET_SHAPE_DERIVATIVES_NATURAL(type) \ - ElementClass::computeDNDS(natural_coord, dnds) - AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); -#undef GET_SHAPE_DERIVATIVES_NATURAL - + auto && dnds = ElementClassHelper<_ek_regular>::getDNDS(natural_coord, type); tangents.mul(dnds, nodes_coord); auto temp_tangents = tangents.transpose(); for (UInt i = 0; i < spatial_dimension - 1; ++i) { auto temp = Vector(temp_tangents(i)); temp_tangents(i) = temp.normalize(); } tangents = temp_tangents.transpose(); } /* -------------------------------------------------------------------------- */ void GeometryUtils::curvature(const Mesh & mesh, const Array & positions, const Element & element, const Vector & natural_coord, Matrix & curvature) { UInt spatial_dimension = mesh.getSpatialDimension(); - auto surface_dimension = spatial_dimension - 1; - const ElementType & type = element.type; + UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); UInt * elem_val = mesh.getConnectivity(type, _not_ghost).storage(); - Matrix dn2ds2(surface_dimension * surface_dimension, - Mesh::getNbNodesPerElement(type)); - -#define GET_SHAPE_SECOND_DERIVATIVES_NATURAL(type) \ - ElementClass::computeDN2DS2(natural_coord, dn2ds2) - AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_SECOND_DERIVATIVES_NATURAL); -#undef GET_SHAPE_SECOND_DERIVATIVES_NATURAL + auto && d2nds2 = + ElementClassHelper<_ek_regular>::getD2NDS2(natural_coord, type); Matrix coords(spatial_dimension, nb_nodes_per_element); mesh.extractNodalValuesFromElement(positions, coords.storage(), elem_val + element.element * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); - curvature.mul(coords, dn2ds2); + curvature.mul(coords, d2nds2); } /* -------------------------------------------------------------------------- */ UInt GeometryUtils::orthogonalProjection( const Mesh & mesh, const Array & positions, const Vector & slave, const Array & elements, Real & gap, Vector & natural_projection, Vector & normal, Real alpha, UInt max_iterations, Real projection_tolerance, Real extension_tolerance) { UInt index = UInt(-1); Real min_gap = std::numeric_limits::max(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt surface_dimension = spatial_dimension - 1; UInt nb_same_sides{0}; UInt nb_boundary_elements{0}; UInt counter{0}; auto & contact_group = mesh.getElementGroup("contact_surface"); for (auto & element : elements) { // filter out elements which are not there in the element group // contact surface created by the surface selector and is stored // in the mesh or mesh_facet, if a element is not there it // returnas UInt(-1) auto & elements_of_type = contact_group.getElements(element.type); if (elements_of_type.find(element.element) == UInt(-1)) continue; // if (!GeometryUtils::isBoundaryElement(mesh, element)) // continue; nb_boundary_elements++; // find the natural coordinate corresponding to the minimum gap // between slave node and master element /* Vector normal_ele(spatial_dimension); GeometryUtils::normal(mesh, positions, element, normal_ele); Vector master(spatial_dimension); GeometryUtils::realProjection(mesh, positions, slave, element, normal_ele, master); Vector xi(natural_projection.size()); GeometryUtils::naturalProjection(mesh, positions, element, master, xi, max_iterations, projection_tolerance);*/ Vector master(spatial_dimension); Vector xi(natural_projection.size()); GeometryUtils::naturalProjection(mesh, positions, element, slave, master, xi, max_iterations, projection_tolerance); Matrix tangent_ele(surface_dimension, spatial_dimension); GeometryUtils::covariantBasis(mesh, positions, element, xi, tangent_ele); Vector normal_ele(spatial_dimension); GeometryUtils::normal(mesh, element, tangent_ele, normal_ele); // if gap between master projection and slave point is zero, then // it means that slave point lies on the master element, hence the // normal from master to slave cannot be computed in that case auto master_to_slave = slave - master; Real temp_gap = master_to_slave.norm(); if (temp_gap != 0) master_to_slave /= temp_gap; // also the slave point should lie inside the master surface in // case of explicit or outside in case of implicit, one way to // check that is by checking the dot product of normal at each // master element, if the values of all dot product is same then // the slave point lies on the same side of each master element // A alpha parameter is introduced which is 1 in case of explicit // and -1 in case of implicit, therefor the variation (dot product // + alpha) should be close to zero (within tolerance) for both // cases Real direction_tolerance = 1e-8; auto product = master_to_slave.dot(normal_ele); auto variation = std::abs(product + alpha); if (variation <= direction_tolerance and temp_gap <= min_gap and GeometryUtils::isValidProjection(xi, extension_tolerance)) { gap = -temp_gap; min_gap = temp_gap; index = counter; natural_projection = xi; normal = normal_ele; } if (temp_gap == 0 or variation <= direction_tolerance) nb_same_sides++; counter++; } // if point is not on the same side of all the boundary elements // than it is not consider even if the closet master element is // found if (nb_same_sides != nb_boundary_elements) index = UInt(-1); return index; } /* -------------------------------------------------------------------------- */ UInt GeometryUtils::orthogonalProjection( const Mesh & mesh, const Array & positions, const Vector & slave, const Array & elements, Real & gap, Vector & natural_projection, Vector & normal, Matrix & tangent, Real /*alpha*/, UInt max_iterations, Real projection_tolerance, Real extension_tolerance) { UInt index = UInt(-1); Real min_gap = std::numeric_limits::max(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt surface_dimension = spatial_dimension - 1; auto & contact_group = mesh.getElementGroup("contact_surface"); for (auto && tuple : enumerate(elements)) { auto & counter = std::get<0>(tuple); auto & element = std::get<1>(tuple); // filter out elements which are not there in the element group // contact surface created by the surface selector and is stored // in the mesh or mesh_facet, if a element is not there it // returnas UInt(-1) auto & elements_of_type = contact_group.getElements(element.type); if (elements_of_type.find(element.element) == UInt(-1)) continue; Vector master(spatial_dimension); Vector xi_ele(natural_projection.size()); GeometryUtils::naturalProjection(mesh, positions, element, slave, master, xi_ele, max_iterations, projection_tolerance); Matrix tangent_ele(surface_dimension, spatial_dimension); GeometryUtils::covariantBasis(mesh, positions, element, xi_ele, tangent_ele); Vector normal_ele(spatial_dimension); GeometryUtils::normal(mesh, element, tangent_ele, normal_ele); // if gap between master projection and slave point is zero, then // it means that slave point lies on the master element, hence the // normal from master to slave cannot be computed in that case auto master_to_slave = slave - master; Real temp_gap = master_to_slave.norm(); if (temp_gap != 0) master_to_slave /= temp_gap; // A alpha parameter is introduced which is 1 in case of explicit // and -1 in case of implicit, therefor the variation (dot product // + alpha) should be close to zero (within tolerance) for both // cases auto product = master_to_slave.dot(normal_ele); if (product < 0 and temp_gap <= min_gap and GeometryUtils::isValidProjection(xi_ele, extension_tolerance)) { gap = -temp_gap; min_gap = temp_gap; index = counter; natural_projection = xi_ele; normal = normal_ele; tangent = tangent_ele; } } return index; } /* -------------------------------------------------------------------------- */ void GeometryUtils::realProjection(const Mesh & mesh, const Array & positions, const Vector & slave, const Element & element, const Vector & normal, Vector & projection) { UInt spatial_dimension = mesh.getSpatialDimension(); const ElementType & type = element.type; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); UInt * elem_val = mesh.getConnectivity(type, _not_ghost).storage(); Matrix nodes_coord(spatial_dimension, nb_nodes_per_element); mesh.extractNodalValuesFromElement(positions, nodes_coord.storage(), elem_val + element.element * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); Vector point(nodes_coord(0)); Real alpha = (slave - point).dot(normal); projection = slave - alpha * normal; } /* -------------------------------------------------------------------------- */ void GeometryUtils::realProjection(const Mesh & mesh, const Array & positions, const Element & element, const Vector & natural_coord, Vector & projection) { UInt spatial_dimension = mesh.getSpatialDimension(); const ElementType & type = element.type; - UInt nb_nodes_per_element = mesh.getNbNodesPerElement(element.type); - - Vector shapes(nb_nodes_per_element); -#define GET_SHAPE_NATURAL(type) \ - ElementClass::computeShapes(natural_coord, shapes) - AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_NATURAL); -#undef GET_SHAPE_NATURAL + auto nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); + auto shapes = + ElementClassHelper<_ek_regular>::getN(natural_coord, element.type); Matrix nodes_coord(spatial_dimension, nb_nodes_per_element); UInt * elem_val = mesh.getConnectivity(type, _not_ghost).storage(); mesh.extractNodalValuesFromElement(positions, nodes_coord.storage(), elem_val + element.element * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); projection.mul(nodes_coord, shapes); } /* -------------------------------------------------------------------------- */ void GeometryUtils::naturalProjection( const Mesh & mesh, const Array & positions, const Element & element, const Vector & slave_coords, Vector & master_coords, Vector & natural_projection, UInt max_iterations, Real projection_tolerance) { UInt spatial_dimension = mesh.getSpatialDimension(); UInt surface_dimension = spatial_dimension - 1; const ElementType & type = element.type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); UInt * elem_val = mesh.getConnectivity(type, _not_ghost).storage(); Matrix nodes_coord(spatial_dimension, nb_nodes_per_element); mesh.extractNodalValuesFromElement(positions, nodes_coord.storage(), elem_val + element.element * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); // initial guess natural_projection.zero(); // obhjective function computed on the natural_guess Matrix f(surface_dimension, 1); // jacobian matrix computed on the natural_guess Matrix J(surface_dimension, surface_dimension); // Jinv = J^{-1} Matrix Jinv(surface_dimension, surface_dimension); // dxi = \xi_{k+1} - \xi_{k} in the iterative process Matrix dxi(surface_dimension, 1); // gradient at natural projection Matrix gradient(surface_dimension, spatial_dimension); // second derivative at natural peojection Matrix double_gradient(surface_dimension, surface_dimension); // second derivative of shape function at natural projection - Matrix dn2ds2(surface_dimension * surface_dimension, + Matrix d2nds2(surface_dimension * surface_dimension, nb_nodes_per_element); auto compute_double_gradient = - [&double_gradient, &dn2ds2, &nodes_coord, surface_dimension, + [&double_gradient, &d2nds2, &nodes_coord, surface_dimension, spatial_dimension](UInt & alpha, UInt & beta) { auto index = alpha * surface_dimension + beta; Vector d_alpha_beta(spatial_dimension); - auto dn2ds2_transpose = dn2ds2.transpose(); - Vector dn2ds2_alpha_beta(dn2ds2_transpose(index)); + auto d2nds2_transpose = d2nds2.transpose(); + Vector d2nds2_alpha_beta(d2nds2_transpose(index)); - d_alpha_beta.mul(nodes_coord, dn2ds2_alpha_beta); + d_alpha_beta.mul(nodes_coord, d2nds2_alpha_beta); return d_alpha_beta; }; /* --------------------------- */ /* init before iteration loop */ /* --------------------------- */ // do interpolation auto update_f = [&f, &master_coords, &natural_projection, &nodes_coord, &slave_coords, &gradient, surface_dimension, spatial_dimension, nb_nodes_per_element, type]() { // compute real coords on natural projection - Vector shapes(nb_nodes_per_element); -#define GET_SHAPE_NATURAL(type) \ - ElementClass::computeShapes(natural_projection, shapes) - AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_NATURAL); -#undef GET_SHAPE_NATURAL + auto && shapes = + ElementClassHelper<_ek_regular>::getN(natural_projection, type); master_coords.mul(nodes_coord, shapes); auto distance = slave_coords - master_coords; // first derivative of shape function at natural projection - Matrix dnds(surface_dimension, nb_nodes_per_element); - - // computing gradient at projection point -#define GET_SHAPE_DERIVATIVES_NATURAL(type) \ - ElementClass::computeDNDS(natural_projection, dnds) - AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); -#undef GET_SHAPE_DERIVATIVES_NATURAL - + auto && dnds = + ElementClassHelper<_ek_regular>::getDNDS(natural_projection, type); gradient.mul(dnds, nodes_coord); // gradient transpose at natural projection Matrix gradient_transpose(surface_dimension, spatial_dimension); gradient_transpose = gradient.transpose(); // loop over surface dimensions for (auto alpha : arange(surface_dimension)) { Vector gradient_alpha(gradient_transpose(alpha)); f(alpha, 0) = -2. * gradient_alpha.dot(distance); } // compute initial error auto error = f.norm(); return error; }; auto projection_error = update_f(); /* --------------------------- */ /* iteration loop */ /* --------------------------- */ UInt iterations{0}; while (projection_tolerance < projection_error and iterations < max_iterations) { // compute covariant components of metric tensor auto a = GeometryUtils::covariantMetricTensor(gradient); // computing second derivative at natural projection -#define GET_SHAPE_SECOND_DERIVATIVES_NATURAL(type) \ - ElementClass::computeDN2DS2(natural_projection, dn2ds2) - AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_SECOND_DERIVATIVES_NATURAL); -#undef GET_SHAPE_SECOND_DERIVATIVES_NATURAL + d2nds2 = ElementClassHelper<_ek_regular>::getD2NDS2(natural_projection, type); // real coord - physical guess auto distance = slave_coords - master_coords; // computing Jacobian J for (auto alpha : arange(surface_dimension)) { for (auto beta : arange(surface_dimension)) { auto dgrad_alpha_beta = compute_double_gradient(alpha, beta); J(alpha, beta) = 2. * (a(alpha, beta) - dgrad_alpha_beta.dot(distance)); } } Jinv.inverse(J); // compute increment dxi.mul(Jinv, f, -1.0); // update our guess natural_projection += Vector(dxi(0)); projection_error = update_f(); iterations++; } /* #define GET_NATURAL_COORDINATE(type) \ ElementClass::inverseMap(slave_coords, nodes_coord, \ natural_projection, max_iterations, projection_tolerance) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NATURAL_COORDINATE); #undef GET_NATURAL_COORDINATE*/ } /* -------------------------------------------------------------------------- */ void GeometryUtils::contravariantBasis(const Matrix & covariant, Matrix & contravariant) { auto inv_A = GeometryUtils::contravariantMetricTensor(covariant); contravariant.mul(inv_A, covariant); } /* -------------------------------------------------------------------------- */ Matrix GeometryUtils::covariantMetricTensor(const Matrix & covariant_bases) { Matrix A(covariant_bases.rows(), covariant_bases.rows()); A.mul(covariant_bases, covariant_bases); return A; } /* -------------------------------------------------------------------------- */ Matrix GeometryUtils::contravariantMetricTensor(const Matrix & covariant_bases) { auto A = GeometryUtils::covariantMetricTensor(covariant_bases); auto inv_A = A.inverse(); return inv_A; } /* -------------------------------------------------------------------------- */ Matrix GeometryUtils::covariantCurvatureTensor( const Mesh & mesh, const Array & positions, const Element & element, const Vector & natural_coord, const Vector & normal) { UInt spatial_dimension = mesh.getSpatialDimension(); auto surface_dimension = spatial_dimension - 1; const ElementType & type = element.type; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt * elem_val = mesh.getConnectivity(type, _not_ghost).storage(); - Matrix dn2ds2(surface_dimension * surface_dimension, - nb_nodes_per_element); - -#define GET_SHAPE_SECOND_DERIVATIVES_NATURAL(type) \ - ElementClass::computeDN2DS2(natural_coord, dn2ds2) - AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_SECOND_DERIVATIVES_NATURAL); -#undef GET_SHAPE_SECOND_DERIVATIVES_NATURAL + auto && d2nds2 = + ElementClassHelper<_ek_regular>::getD2NDS2(natural_coord, type); Matrix coords(spatial_dimension, nb_nodes_per_element); mesh.extractNodalValuesFromElement(positions, coords.storage(), elem_val + element.element * nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); Matrix curvature(spatial_dimension, surface_dimension * surface_dimension); - curvature.mul(coords, dn2ds2); + curvature.mul(coords, d2nds2); Matrix H(surface_dimension, surface_dimension); UInt i = 0; for (auto alpha : arange(surface_dimension)) { for (auto beta : arange(surface_dimension)) { Vector temp(curvature(i)); H(alpha, beta) = temp.dot(normal); i++; } } return H; } } // namespace akantu diff --git a/src/model/contact_mechanics/resolution_utils.cc b/src/model/contact_mechanics/resolution_utils.cc index 9c8670571..4530f4901 100644 --- a/src/model/contact_mechanics/resolution_utils.cc +++ b/src/model/contact_mechanics/resolution_utils.cc @@ -1,330 +1,327 @@ /** * @file resolution_utils.cc * * @author Mohit Pundir * * @date creation: Mon Mmay 20 2019 * @date last modification: Mon May 20 2019 * * @brief Implementation of various utilities neede for resolution class * * @section LICENSE * * Copyright (©) 2010-2018 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 "resolution_utils.hh" +#include "element_class_helper.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ void ResolutionUtils::computeShapeFunctionMatric( const ContactElement & element, const Vector & projection, Matrix & shape_matric) { shape_matric.zero(); const ElementType & type = element.master.type; auto surface_dimension = Mesh::getSpatialDimension(type); auto spatial_dimension = surface_dimension + 1; UInt nb_nodes_per_contact = element.getNbNodes(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); AKANTU_DEBUG_ASSERT(spatial_dimension == shape_matric.rows() && spatial_dimension * nb_nodes_per_contact == shape_matric.cols(), "Shape Matric dimensions are not correct"); - Vector shapes(nb_nodes_per_element); - -#define GET_SHAPE_NATURAL(type) \ - ElementClass::computeShapes(projection, shapes) - AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_NATURAL); -#undef GET_SHAPE_NATURAL + auto && shapes = + ElementClassHelper<_ek_regular>::getN(projection, type); for (auto i : arange(nb_nodes_per_contact)) { for (auto j : arange(spatial_dimension)) { if (i == 0) { shape_matric(j, i * spatial_dimension + j) = 1; continue; } shape_matric(j, i * spatial_dimension + j) = -shapes[i - 1]; } } } /* -------------------------------------------------------------------------- */ /*void ResolutionUtils::firstVariationNormalGap(const ContactElement & element, const Vector & projection, const Vector & normal, Vector & delta_g) { delta_g.clear(); const auto & type = element.master.type; auto surface_dimension = Mesh::getSpatialDimension(type); auto spatial_dimension = surface_dimension + 1; Vector shapes(Mesh::getNbNodesPerElement(type)); #define GET_SHAPES_NATURAL(type) \ ElementClass::computeShapes(projection, shapes) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPES_NATURAL); #undef GET_SHAPES_NATURAL for (UInt i : arange(spatial_dimension)) { delta_g[i] = normal[i]; for (UInt j : arange(shapes.size())) { delta_g[(1 + j) * spatial_dimension + i] = -shapes[j] * normal[i]; } } }*/ /* -------------------------------------------------------------------------- */ /*void ResolutionUtils::secondVariationNormalGap(const ContactElement & element, const Matrix & covariant_basis, const Matrix & curvature, const Vector & projection, const Vector & normal, Real & gap, Matrix & ddelta_g) { const auto & type = element.master.type; auto surface_dimension = Mesh::getSpatialDimension(type); auto spatial_dimension = surface_dimension + 1; UInt nb_nodes = element.getNbNodes(); Array dnds_n(nb_nodes * spatial_dimension, surface_dimension); ResolutionUtils::computeNalpha(element, projection, normal, dnds_n); Array delta_xi(nb_nodes * spatial_dimension, surface_dimension); ResolutionUtils::firstVariationNaturalCoordinate(element, covariant_basis, projection, normal, gap, delta_xi); Matrix a_alpha_beta(surface_dimension, surface_dimension); ResolutionUtils::computeMetricTensor(a_alpha_beta, covariant_basis); a_alpha_beta = a_alpha_beta.inverse(); Matrix h_alpha_beta(surface_dimension, surface_dimension); ResolutionUtils::computeSecondMetricTensor(element, curvature, normal, h_alpha_beta); for (auto && values : zip(arange(surface_dimension), make_view(dnds_n, dnds_n.size()), make_view(delta_xi, delta_xi.size()))) { auto & alpha = std::get<0>(values); auto & dnds_n_alpha = std::get<1>(values); auto & delta_xi_alpha = std::get<2>(values); // term 1 from Numerical methods in contact mechanics : Vlad // Yastrebov eq 2.48 Matrix mat_n(dnds_n_alpha.storage(), dnds_n_alpha.size(), 1); Matrix mat_xi(delta_xi_alpha.storage(), delta_xi_alpha.size(), 1); Matrix tmp1(dnds_n_alpha.size(), dnds_n_alpha.size()); tmp1.mul(mat_n, mat_xi, -1); Matrix tmp2(dnds_n_alpha.size(), dnds_n_alpha.size()); tmp2.mul(mat_xi, mat_n, -1); Matrix term1(dnds_n_alpha.size(), dnds_n_alpha.size()); term1 = tmp1 + tmp2; // computing term 2 & term 3 from Numerical methods in contact // mechanics : Vlad Yastrebov eq 2.48 Matrix term2(delta_xi_alpha.size(), delta_xi_alpha.size()); Matrix term3(dnds_n_alpha.size(), dnds_n_alpha.size()); for (auto && values2 : zip(arange(surface_dimension), make_view(dnds_n, dnds_n.size()), make_view(delta_xi, delta_xi.size()))) { auto & beta = std::get<0>(values2); auto & dnds_n_beta = std::get<1>(values2); auto & delta_xi_beta = std::get<2>(values2); // term 2 Matrix mat_xi_beta(delta_xi_beta.storage(), delta_xi.size(), 1); Matrix tmp3(delta_xi_beta.size(), delta_xi_beta.size()); Real pre_factor = h_alpha_beta(alpha, beta); for (auto k : arange(surface_dimension)) { for (auto m : arange(surface_dimension)) { pre_factor -= gap * h_alpha_beta(alpha, k) * a_alpha_beta(k, m) * h_alpha_beta(m, beta); } } pre_factor *= -1.; tmp3.mul(mat_xi, mat_xi_beta, pre_factor); // term 3 Matrix mat_n_beta(dnds_n_beta.storage(), dnds_n_beta.size(), 1); Real factor = gap * a_alpha_beta(alpha, beta); Matrix tmp4(dnds_n_alpha.size(), dnds_n_alpha.size()); tmp4.mul(mat_n, mat_n_beta, factor); term3 += tmp4; } ddelta_g += term1 + term2 + term3; } }*/ /* -------------------------------------------------------------------------- */ /*void ResolutionUtils::computeTalpha(const ContactElement & element, const Matrix & covariant_basis, const Vector & projection, Array & t_alpha) { t_alpha.clear(); const auto & type = element.master.type; auto surface_dimension = Mesh::getSpatialDimension(type); auto spatial_dimension = surface_dimension + 1; Vector shapes(Mesh::getNbNodesPerElement(type)); #define GET_SHAPES_NATURAL(type) \ ElementClass::computeShapes(projection, shapes) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPES_NATURAL); #undef GET_SHAPES_NATURAL for (auto && values : zip(covariant_basis.transpose(), make_view(t_alpha, t_alpha.size()))) { auto & tangent_beta = std::get<0>(values); auto & t_beta = std::get<1>(values); for (UInt i : arange(spatial_dimension)) { t_beta[i] = tangent_beta(i); for (UInt j : arange(shapes.size())) { t_beta[(1 + j) * spatial_dimension + i] = -shapes[j] * tangent_beta(i); } } } }*/ /* -------------------------------------------------------------------------- */ /*void ResolutionUtils::computeNalpha(const ContactElement & element, const Vector & projection, const Vector & normal, Array & n_alpha) { n_alpha.clear(); const auto & type = element.master.type; auto surface_dimension = Mesh::getSpatialDimension(type); auto spatial_dimension = surface_dimension + 1; Matrix shape_derivatives(surface_dimension, Mesh::getNbNodesPerElement(type)); #define GET_SHAPE_DERIVATIVES_NATURAL(type) \ ElementClass::computeDNDS(projection, shape_derivatives) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); #undef GET_SHAPE_DERIVATIVES_NATURAL for (auto && values : zip(shape_derivatives.transpose(), make_view(n_alpha, n_alpha.size()))) { auto & dnds = std::get<0>(values); auto & n_s = std::get<1>(values); for (UInt i : arange(spatial_dimension)) { n_s[i] = 0; for (UInt j : arange(dnds.size())) { n_s[(1 + j) * spatial_dimension + i] = -dnds(j) * normal[i]; } } } }*/ /* -------------------------------------------------------------------------- */ /*void ResolutionUtils::firstVariationNaturalCoordinate(const ContactElement & element, const Matrix & covariant_basis, const Vector & projection, const Vector & normal, const Real & gap, Array & delta_xi) { delta_xi.clear(); const auto & type = element.master.type; auto surface_dimension = Mesh::getSpatialDimension(type); auto spatial_dimension = surface_dimension + 1; auto inv_A = GeometryUtils::contravariantMetricTensor(covariant_basis); auto nb_nodes = element.getNbNodes(); Array t_alpha(nb_nodes * spatial_dimension, surface_dimension); Array n_alpha(nb_nodes * spatial_dimension, surface_dimension); ResolutionUtils::computeTalpha(element, covariant_basis, projection, t_alpha); ResolutionUtils::computeNalpha(element, projection, normal, n_alpha); for (auto && entry : zip(arange(surface_dimension), make_view(delta_xi, delta_xi.size()))) { auto & alpha = std::get<0>(entry); auto & d_alpha = std::get<1>(entry); for (auto && values : zip(arange(surface_dimension), make_view(t_alpha, t_alpha.size()), make_view(n_alpha, n_alpha.size()))) { auto & beta = std::get<0>(values); auto & t_beta = std::get<1>(values); //auto & n_beta = std::get<2>(values); //d_alpha += (t_beta + gap * n_beta) * m_alpha_beta(alpha, //beta); d_alpha += t_beta * inv_A(alpha, beta); } } }*/ /* -------------------------------------------------------------------------- */ /*void ResolutionUtils::computeMetricTensor(Matrix & m_alpha_beta, const Matrix & tangents) { m_alpha_beta.mul(tangents, tangents); }*/ /* -------------------------------------------------------------------------- */ /*void ResolutionUtils::computeSecondMetricTensor(const ContactElement & element, const Matrix & curvature, const Vector & normal, Matrix & metric) { const auto & type = element.master.type; auto surface_dimension = Mesh::getSpatialDimension(type); auto i = 0; for (auto alpha : arange(surface_dimension) ) { for (auto beta : arange(surface_dimension)) { Vector temp(curvature(i)); metric(alpha, beta) = normal.dot(temp); i++; } } }*/ } // namespace akantu diff --git a/src/model/contact_mechanics/resolutions/resolution_penalty.cc b/src/model/contact_mechanics/resolutions/resolution_penalty.cc index 598aff62f..cadd1db32 100644 --- a/src/model/contact_mechanics/resolutions/resolution_penalty.cc +++ b/src/model/contact_mechanics/resolutions/resolution_penalty.cc @@ -1,903 +1,867 @@ /** * @file resolution_penalty.cc * * @author Mohit Pundir * * @date creation: Mon Jan 7 2019 * @date last modification: Mon Jan 7 2019 * * @brief Specialization of the resolution class for the penalty method * * @section LICENSE * * Copyright (©) 2010-2018 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 "resolution_penalty.hh" +#include "element_class_helper.hh" +/* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ ResolutionPenalty::ResolutionPenalty(ContactMechanicsModel & model, const ID & id) : Resolution(model, id) { AKANTU_DEBUG_IN(); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::initialize() { this->registerParam("epsilon_n", epsilon_n, Real(0.), _pat_parsable | _pat_modifiable, "Normal penalty parameter"); this->registerParam("epsilon_t", epsilon_t, Real(0.), _pat_parsable | _pat_modifiable, "Tangential penalty parameter"); } /* -------------------------------------------------------------------------- */ Real ResolutionPenalty::computeNormalTraction(Real & gap) { return epsilon_n * macaulay(gap); } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeNormalForce(const ContactElement & element, Vector & force) { force.zero(); auto & gaps = model.getGaps(); auto & projections = model.getProjections(); auto & normals = model.getNormals(); auto surface_dimension = spatial_dimension - 1; Real gap(gaps.begin()[element.slave]); Vector normal(normals.begin(spatial_dimension)[element.slave]); Vector projection(projections.begin(surface_dimension)[element.slave]); auto & nodal_area = const_cast &>(model.getNodalArea()); // compute normal traction Real p_n = computeNormalTraction(gap); p_n *= nodal_area[element.slave]; UInt nb_nodes_per_contact = element.getNbNodes(); Matrix shape_matric(spatial_dimension, spatial_dimension * nb_nodes_per_contact); ResolutionUtils::computeShapeFunctionMatric(element, projection, shape_matric); force.mul(shape_matric, normal, p_n); } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeTangentialForce(const ContactElement & element, Vector & force) { if (mu == 0) return; force.zero(); UInt surface_dimension = spatial_dimension - 1; // compute covariant basis auto & projections = model.getProjections(); Vector projection(projections.begin(surface_dimension)[element.slave]); auto & normals = model.getNormals(); Vector normal(normals.begin(spatial_dimension)[element.slave]); auto & tangents = model.getTangents(); Matrix covariant_basis( tangents.begin(surface_dimension, spatial_dimension)[element.slave]); // check for no-contact to contact condition // need a better way to check if new node added is not presnt in the // previous master elemets auto & previous_master_elements = model.getPreviousMasterElements(); if (element.slave >= previous_master_elements.size()) return; auto & previous_element = previous_master_elements[element.slave]; if (previous_element.type == _not_defined) return; // compute tangential traction using return map algorithm auto & tangential_tractions = model.getTangentialTractions(); Vector tangential_traction( tangential_tractions.begin(surface_dimension)[element.slave]); this->computeTangentialTraction(element, covariant_basis, tangential_traction); UInt nb_nodes_per_contact = element.getNbNodes(); Matrix shape_matric(spatial_dimension, spatial_dimension * nb_nodes_per_contact); ResolutionUtils::computeShapeFunctionMatric(element, projection, shape_matric); auto contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(covariant_basis); auto & nodal_area = const_cast &>(model.getNodalArea()); for (auto && values1 : enumerate(covariant_basis.transpose())) { auto & alpha = std::get<0>(values1); auto & tangent_alpha = std::get<1>(values1); for (auto && values2 : enumerate(tangential_traction)) { auto & beta = std::get<0>(values2); auto & traction_beta = std::get<1>(values2); Vector tmp(force.size()); tmp.mul(shape_matric, tangent_alpha, traction_beta); tmp *= contravariant_metric_tensor(alpha, beta) * nodal_area[element.slave]; force += tmp; } } } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeTangentialTraction( const ContactElement & element, const Matrix & covariant_basis, Vector & traction_tangential) { UInt surface_dimension = spatial_dimension - 1; auto & gaps = model.getGaps(); auto & gap = gaps.begin()[element.slave]; // Return map algorithm is employed // compute trial traction Vector traction_trial(surface_dimension); this->computeTrialTangentialTraction(element, covariant_basis, traction_trial); // compute norm of trial traction Real traction_trial_norm = 0; auto contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(covariant_basis); for (auto i : arange(surface_dimension)) { for (auto j : arange(surface_dimension)) { traction_trial_norm += traction_trial[i] * traction_trial[j] * contravariant_metric_tensor(i, j); } } traction_trial_norm = sqrt(traction_trial_norm); // check stick or slip condition auto & contact_state = model.getContactState(); auto & state = contact_state.begin()[element.slave]; Real p_n = computeNormalTraction(gap); bool stick = (traction_trial_norm <= mu * p_n) ? true : false; if (stick) { state = ContactState::_stick; computeStickTangentialTraction(element, traction_trial, traction_tangential); } else { state = ContactState::_slip; computeSlipTangentialTraction(element, covariant_basis, traction_trial, traction_tangential); } } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeTrialTangentialTraction( const ContactElement & element, const Matrix & covariant_basis, Vector & traction) { UInt surface_dimension = spatial_dimension - 1; auto & projections = model.getProjections(); Vector current_projection( projections.begin(surface_dimension)[element.slave]); auto & previous_projections = model.getPreviousProjections(); Vector previous_projection( previous_projections.begin(surface_dimension)[element.slave]); // method from Laursen et. al. /*auto covariant_metric_tensor = GeometryUtils::covariantMetricTensor(covariant_basis); auto increment_projection = current_projection - previous_projection; traction.mul(covariant_metric_tensor, increment_projection, epsilon_t); auto & previous_tangential_tractions = model.getPreviousTangentialTractions(); Vector previous_traction(previous_tangential_tractions.begin(surface_dimension)[element.slave]); traction = previous_traction + traction;*/ // method from Schweizerhof auto covariant_metric_tensor = GeometryUtils::covariantMetricTensor(covariant_basis); auto & previous_tangential_tractions = model.getPreviousTangentialTractions(); Vector previous_traction( previous_tangential_tractions.begin(surface_dimension)[element.slave]); auto & previous_tangents = model.getPreviousTangents(); Matrix previous_covariant_basis(previous_tangents.begin( surface_dimension, spatial_dimension)[element.slave]); auto previous_contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(previous_covariant_basis); auto current_tangent = covariant_basis.transpose(); auto previous_tangent = previous_covariant_basis.transpose(); for (auto alpha : arange(surface_dimension)) { Vector tangent_alpha(current_tangent(alpha)); for (auto gamma : arange(surface_dimension)) { for (auto beta : arange(surface_dimension)) { Vector tangent_beta(previous_tangent(beta)); auto t_alpha_t_beta = tangent_beta.dot(tangent_alpha); traction[alpha] += previous_traction[gamma] * previous_contravariant_metric_tensor(gamma, beta) * t_alpha_t_beta; } } } auto & previous_master_elements = model.getPreviousMasterElements(); auto & previous_element = previous_master_elements[element.slave]; Vector previous_real_projection(spatial_dimension); GeometryUtils::realProjection( model.getMesh(), model.getContactDetector().getPositions(), previous_element, previous_projection, previous_real_projection); Vector current_real_projection(spatial_dimension); GeometryUtils::realProjection( model.getMesh(), model.getContactDetector().getPositions(), element.master, current_projection, current_real_projection); auto increment_real = current_real_projection - previous_real_projection; Vector increment_xi(surface_dimension); auto contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(covariant_basis); // increment in natural coordinate for (auto beta : arange(surface_dimension)) { for (auto gamma : arange(surface_dimension)) { auto temp = increment_real.dot(current_tangent(gamma)); temp *= contravariant_metric_tensor(beta, gamma); increment_xi[beta] += temp; } } Vector temp(surface_dimension); temp.mul(covariant_metric_tensor, increment_xi, epsilon_t); traction -= temp; } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeStickTangentialTraction( const ContactElement & /*element*/, Vector & traction_trial, Vector & traction_tangential) { traction_tangential = traction_trial; } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeSlipTangentialTraction( const ContactElement & element, const Matrix & covariant_basis, Vector & traction_trial, Vector & traction_tangential) { UInt surface_dimension = spatial_dimension - 1; auto & gaps = model.getGaps(); auto & gap = gaps.begin()[element.slave]; // compute norm of trial traction Real traction_trial_norm = 0; auto contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(covariant_basis); for (auto alpha : arange(surface_dimension)) { for (auto beta : arange(surface_dimension)) { traction_trial_norm += traction_trial[alpha] * traction_trial[beta] * contravariant_metric_tensor(alpha, beta); } } traction_trial_norm = sqrt(traction_trial_norm); auto slip_direction = traction_trial; slip_direction /= traction_trial_norm; Real p_n = computeNormalTraction(gap); traction_tangential = slip_direction; traction_tangential *= mu * p_n; } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeNormalModuli(const ContactElement & element, Matrix & stiffness) { auto surface_dimension = spatial_dimension - 1; auto & gaps = model.getGaps(); Real gap(gaps.begin()[element.slave]); auto & projections = model.getProjections(); Vector projection(projections.begin(surface_dimension)[element.slave]); auto & nodal_areas = model.getNodalArea(); auto & nodal_area = nodal_areas.begin()[element.slave]; auto & normals = model.getNormals(); Vector normal(normals.begin(spatial_dimension)[element.slave]); - auto & mesh = model.getMesh(); - // method from Schweizerhof and A. Konyukhov, K. Schweizerhof // DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3 // construct A matrix const ElementType & type = element.master.type; - UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); - Vector shapes(nb_nodes_per_element); -#define GET_SHAPE_NATURAL(type) \ - ElementClass::computeShapes(projection, shapes) - AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_NATURAL); -#undef GET_SHAPE_NATURAL + auto && shapes = ElementClassHelper<_ek_regular>::getN(projection, type); UInt nb_nodes_per_contact = element.getNbNodes(); Matrix A(spatial_dimension, spatial_dimension * nb_nodes_per_contact); for (auto i : arange(nb_nodes_per_contact)) { for (auto j : arange(spatial_dimension)) { if (i == 0) { A(j, i * spatial_dimension + j) = 1; continue; } A(j, i * spatial_dimension + j) = -shapes[i - 1]; } } // construct the main part of normal matrix Matrix k_main(nb_nodes_per_contact * spatial_dimension, nb_nodes_per_contact * spatial_dimension); Matrix n_outer_n(spatial_dimension, spatial_dimension); Matrix mat_n(normal.storage(), normal.size(), 1.); n_outer_n.mul(mat_n, mat_n); Matrix tmp(spatial_dimension, spatial_dimension * nb_nodes_per_contact); tmp.mul(n_outer_n, A); k_main.mul(A, tmp); k_main *= epsilon_n * heaviside(gap) * nodal_area; // construct the rotational part of the normal matrix auto & tangents = model.getTangents(); Matrix covariant_basis( tangents.begin(surface_dimension, spatial_dimension)[element.slave]); // GeometryUtils::covariantBasis(model.getMesh(), // model.getContactDetector().getPositions(), element.master, // normal, // projection, covariant_basis); auto contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(covariant_basis); // computing shape derivatives - Matrix shape_derivatives(surface_dimension, - Mesh::getNbNodesPerElement(type)); - -#define GET_SHAPE_DERIVATIVES_NATURAL(type) \ - ElementClass::computeDNDS(projection, shape_derivatives) - AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); -#undef GET_SHAPE_DERIVATIVES_NATURAL + auto && shape_derivatives = + ElementClassHelper<_ek_regular>::getDNDS(projection, type); // consists of 2 rotational parts Matrix k_rot1(nb_nodes_per_contact * spatial_dimension, nb_nodes_per_contact * spatial_dimension); Matrix k_rot2(nb_nodes_per_contact * spatial_dimension, nb_nodes_per_contact * spatial_dimension); Matrix Aj(spatial_dimension, spatial_dimension * nb_nodes_per_contact); auto construct_Aj = [&](auto && dnds) { for (auto i : arange(nb_nodes_per_contact)) { for (auto j : arange(spatial_dimension)) { if (i == 0) { Aj(j, i * spatial_dimension + j) = 0; continue; } Aj(j, i * spatial_dimension + j) = dnds(i - 1); } } }; for (auto && values1 : enumerate(covariant_basis.transpose())) { auto & alpha = std::get<0>(values1); auto & tangent = std::get<1>(values1); Matrix n_outer_t(spatial_dimension, spatial_dimension); Matrix mat_t(tangent.storage(), tangent.size(), 1.); n_outer_t.mul(mat_n, mat_t); Matrix t_outer_n(spatial_dimension, spatial_dimension); t_outer_n.mul(mat_t, mat_n); for (auto && values2 : enumerate(shape_derivatives.transpose())) { auto & beta = std::get<0>(values2); auto & dnds = std::get<1>(values2); // construct Aj from shape function wrt to jth natural // coordinate construct_Aj(dnds); Matrix tmp(spatial_dimension, spatial_dimension * nb_nodes_per_contact); Matrix tmp1(nb_nodes_per_contact * spatial_dimension, spatial_dimension * nb_nodes_per_contact); tmp.mul(n_outer_t, A); tmp1.mul(Aj, tmp); tmp1 *= contravariant_metric_tensor(alpha, beta); k_rot1 += tmp1; tmp.mul(t_outer_n, Aj); tmp1.mul(A, tmp); tmp1 *= contravariant_metric_tensor(alpha, beta); k_rot2 += tmp1; } } k_rot1 *= -epsilon_n * heaviside(gap) * gap * nodal_area; k_rot2 *= -epsilon_n * heaviside(gap) * gap * nodal_area; stiffness += k_main + k_rot1 + k_rot2; } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeTangentialModuli(const ContactElement & element, Matrix & stiffness) { if (mu == 0) return; stiffness.zero(); auto & contact_state = model.getContactState(); UInt state = contact_state.begin()[element.slave]; switch (state) { case ContactState::_stick: { computeStickModuli(element, stiffness); break; } case ContactState::_slip: { computeSlipModuli(element, stiffness); break; } default: break; } } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeStickModuli(const ContactElement & element, Matrix & stiffness) { auto surface_dimension = spatial_dimension - 1; auto & projections = model.getProjections(); Vector projection(projections.begin(surface_dimension)[element.slave]); auto & nodal_areas = model.getNodalArea(); auto & nodal_area = nodal_areas.begin()[element.slave]; - auto & mesh = model.getMesh(); - // method from Schweizerhof and A. Konyukhov, K. Schweizerhof // DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3 // construct A matrix const ElementType & type = element.master.type; - UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); - Vector shapes(nb_nodes_per_element); - -#define GET_SHAPE_NATURAL(type) \ - ElementClass::computeShapes(projection, shapes) - AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_NATURAL); -#undef GET_SHAPE_NATURAL + auto && shapes = ElementClassHelper<_ek_regular>::getN(projection, type); UInt nb_nodes_per_contact = element.getNbNodes(); Matrix A(spatial_dimension, spatial_dimension * nb_nodes_per_contact); for (auto i : arange(nb_nodes_per_contact)) { for (auto j : arange(spatial_dimension)) { if (i == 0) { A(j, i * spatial_dimension + j) = 1; continue; } A(j, i * spatial_dimension + j) = -shapes[i - 1]; } } // computing shape derivatives - Matrix shape_derivatives(surface_dimension, - Mesh::getNbNodesPerElement(type)); - -#define GET_SHAPE_DERIVATIVES_NATURAL(type) \ - ElementClass::computeDNDS(projection, shape_derivatives) - AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); -#undef GET_SHAPE_DERIVATIVES_NATURAL + auto && shape_derivatives = + ElementClassHelper<_ek_regular>::getDNDS(projection, type); Matrix Aj(spatial_dimension, spatial_dimension * nb_nodes_per_contact); auto construct_Aj = [&](auto && dnds) { for (auto i : arange(nb_nodes_per_contact)) { for (auto j : arange(spatial_dimension)) { if (i == 0) { Aj(j, i * spatial_dimension + j) = 0; continue; } Aj(j, i * spatial_dimension + j) = dnds(i - 1); } } }; // tangents should have been calculated in normal modulii auto & tangents = model.getTangents(); Matrix covariant_basis( tangents.begin(surface_dimension, spatial_dimension)[element.slave]); auto contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(covariant_basis); // construct 1st part of the stick modulii Matrix k_main(nb_nodes_per_contact * spatial_dimension, nb_nodes_per_contact * spatial_dimension); for (auto && values1 : enumerate(covariant_basis.transpose())) { auto & alpha = std::get<0>(values1); auto & tangent_alpha = std::get<1>(values1); Matrix t_outer_t(spatial_dimension, spatial_dimension); Matrix mat_t_alpha(tangent_alpha.storage(), tangent_alpha.size(), 1.); for (auto && values2 : enumerate(covariant_basis.transpose())) { auto & beta = std::get<0>(values2); auto & tangent_beta = std::get<1>(values2); Matrix mat_t_beta(tangent_beta.storage(), tangent_beta.size(), 1.); t_outer_t.mul(mat_t_alpha, mat_t_beta); Matrix tmp(spatial_dimension, spatial_dimension * nb_nodes_per_contact); Matrix tmp1(nb_nodes_per_contact * spatial_dimension, spatial_dimension * nb_nodes_per_contact); tmp.mul(t_outer_t, A); tmp1.mul(A, tmp); tmp1 *= contravariant_metric_tensor(alpha, beta); k_main += tmp1; } } k_main *= -epsilon_t; // construct 2nd part of the stick modulii auto & tangential_tractions = model.getTangentialTractions(); Vector tangential_traction( tangential_tractions.begin(surface_dimension)[element.slave]); Matrix k_second(nb_nodes_per_contact * spatial_dimension, nb_nodes_per_contact * spatial_dimension); for (auto alpha : arange(surface_dimension)) { Matrix k_sum(nb_nodes_per_contact * spatial_dimension, nb_nodes_per_contact * spatial_dimension); for (auto && values1 : enumerate(shape_derivatives.transpose())) { auto & beta = std::get<0>(values1); auto & dnds = std::get<1>(values1); // construct Aj from shape function wrt to jth natural // coordinate construct_Aj(dnds); for (auto && values2 : enumerate(covariant_basis.transpose())) { auto & gamma = std::get<0>(values2); auto & tangent_gamma = std::get<1>(values2); Matrix t_outer_t(spatial_dimension, spatial_dimension); Matrix mat_t_gamma(tangent_gamma.storage(), tangent_gamma.size(), 1.); for (auto && values3 : enumerate(covariant_basis.transpose())) { auto & theta = std::get<0>(values3); auto & tangent_theta = std::get<1>(values3); Matrix mat_t_theta(tangent_theta.storage(), tangent_theta.size(), 1.); t_outer_t.mul(mat_t_gamma, mat_t_theta); Matrix tmp(spatial_dimension, spatial_dimension * nb_nodes_per_contact); Matrix tmp1(nb_nodes_per_contact * spatial_dimension, spatial_dimension * nb_nodes_per_contact); tmp.mul(t_outer_t, Aj); tmp1.mul(A, tmp); tmp1 *= contravariant_metric_tensor(alpha, theta) * contravariant_metric_tensor(beta, gamma); Matrix tmp2(spatial_dimension, spatial_dimension * nb_nodes_per_contact); Matrix tmp3(nb_nodes_per_contact * spatial_dimension, spatial_dimension * nb_nodes_per_contact); tmp2.mul(t_outer_t, A); tmp3.mul(Aj, tmp2); tmp3 *= contravariant_metric_tensor(alpha, gamma) * contravariant_metric_tensor(beta, theta); k_sum += tmp1 + tmp3; } } } k_second += tangential_traction[alpha] * k_sum; } stiffness += k_main * nodal_area - k_second * nodal_area; } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::computeSlipModuli(const ContactElement & element, Matrix & stiffness) { auto surface_dimension = spatial_dimension - 1; auto & gaps = model.getGaps(); Real gap(gaps.begin()[element.slave]); auto & nodal_areas = model.getNodalArea(); auto & nodal_area = nodal_areas.begin()[element.slave]; // compute normal traction Real p_n = computeNormalTraction(gap); auto & projections = model.getProjections(); Vector projection(projections.begin(surface_dimension)[element.slave]); auto & normals = model.getNormals(); Vector normal(normals.begin(spatial_dimension)[element.slave]); // restructure normal as a matrix for an outer product Matrix mat_n(normal.storage(), normal.size(), 1.); - auto & mesh = model.getMesh(); - // method from Schweizerhof and A. Konyukhov, K. Schweizerhof // DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3 // construct A matrix const ElementType & type = element.master.type; - UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); - Vector shapes(nb_nodes_per_element); - -#define GET_SHAPE_NATURAL(type) \ - ElementClass::computeShapes(projection, shapes) - AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_NATURAL); -#undef GET_SHAPE_NATURAL + auto && shapes = ElementClassHelper<_ek_regular>::getN(projection, type); UInt nb_nodes_per_contact = element.getNbNodes(); Matrix A(spatial_dimension, spatial_dimension * nb_nodes_per_contact); for (auto i : arange(nb_nodes_per_contact)) { for (auto j : arange(spatial_dimension)) { if (i == 0) { A(j, i * spatial_dimension + j) = 1; continue; } A(j, i * spatial_dimension + j) = -shapes[i - 1]; } } // computing shape derivatives - Matrix shape_derivatives(surface_dimension, - Mesh::getNbNodesPerElement(type)); - -#define GET_SHAPE_DERIVATIVES_NATURAL(type) \ - ElementClass::computeDNDS(projection, shape_derivatives) - AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); -#undef GET_SHAPE_DERIVATIVES_NATURAL + auto && shape_derivatives = + ElementClassHelper<_ek_regular>::getDNDS(projection, type); Matrix Aj(spatial_dimension, spatial_dimension * nb_nodes_per_contact); auto construct_Aj = [&](auto && dnds) { for (auto i : arange(nb_nodes_per_contact)) { for (auto j : arange(spatial_dimension)) { if (i == 0) { Aj(j, i * spatial_dimension + j) = 0; continue; } Aj(j, i * spatial_dimension + j) = dnds(i - 1); } } }; // tangents should have been calculated in normal modulii auto & tangents = model.getTangents(); Matrix covariant_basis( tangents.begin(surface_dimension, spatial_dimension)[element.slave]); auto & tangential_tractions = model.getTangentialTractions(); Vector tangential_traction( tangential_tractions.begin(surface_dimension)[element.slave]); // compute norm of trial traction Real traction_norm = 0; auto contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(covariant_basis); for (auto i : arange(surface_dimension)) { for (auto j : arange(surface_dimension)) { traction_norm += tangential_traction[i] * tangential_traction[j] * contravariant_metric_tensor(i, j); } } traction_norm = sqrt(traction_norm); // construct four parts of stick modulii (eq 107,107a-c) Matrix k_first(nb_nodes_per_contact * spatial_dimension, nb_nodes_per_contact * spatial_dimension); Matrix k_second(nb_nodes_per_contact * spatial_dimension, nb_nodes_per_contact * spatial_dimension); Matrix k_third(nb_nodes_per_contact * spatial_dimension, nb_nodes_per_contact * spatial_dimension); Matrix k_fourth(nb_nodes_per_contact * spatial_dimension, nb_nodes_per_contact * spatial_dimension); for (auto && values1 : enumerate(covariant_basis.transpose())) { auto & alpha = std::get<0>(values1); auto & tangent_alpha = std::get<1>(values1); Matrix mat_t_alpha(tangent_alpha.storage(), tangent_alpha.size(), 1.); Matrix t_outer_n(spatial_dimension, spatial_dimension); Matrix t_outer_t(spatial_dimension, spatial_dimension); for (auto && values2 : zip(arange(surface_dimension), covariant_basis.transpose(), shape_derivatives.transpose())) { auto & beta = std::get<0>(values2); auto & tangent_beta = std::get<1>(values2); auto & dnds = std::get<2>(values2); // construct Aj from shape function wrt to jth natural // coordinate construct_Aj(dnds); // eq 107 Matrix mat_t_beta(tangent_beta.storage(), tangent_beta.size(), 1.); t_outer_n.mul(mat_t_beta, mat_n); Matrix tmp(spatial_dimension, spatial_dimension * nb_nodes_per_contact); Matrix tmp1(nb_nodes_per_contact * spatial_dimension, spatial_dimension * nb_nodes_per_contact); tmp.mul(t_outer_n, A); tmp1.mul(A, tmp); tmp1 *= epsilon_n * mu * tangential_traction[alpha] * contravariant_metric_tensor(alpha, beta); tmp1 /= traction_norm; k_first += tmp1 * nodal_area; // eq 107a t_outer_t.mul(mat_t_alpha, mat_t_beta); tmp.mul(t_outer_t, A); tmp1.mul(A, tmp); tmp1 *= epsilon_t * mu * p_n * contravariant_metric_tensor(alpha, beta); tmp1 /= traction_norm; k_second += tmp1 * nodal_area; for (auto && values3 : enumerate(covariant_basis.transpose())) { auto & gamma = std::get<0>(values3); auto & tangent_gamma = std::get<1>(values3); Matrix mat_t_gamma(tangent_gamma.storage(), tangent_gamma.size(), 1.); for (auto && values4 : enumerate(covariant_basis.transpose())) { auto & theta = std::get<0>(values4); auto & tangent_theta = std::get<1>(values4); Matrix mat_t_theta(tangent_theta.storage(), tangent_theta.size(), 1.); t_outer_t.mul(mat_t_gamma, mat_t_theta); // eq 107b tmp.mul(t_outer_t, A); tmp1.mul(A, tmp); tmp1 *= epsilon_t * mu * p_n * tangential_traction[alpha] * tangential_traction[beta]; tmp1 *= contravariant_metric_tensor(alpha, gamma) * contravariant_metric_tensor(beta, theta); tmp1 /= pow(traction_norm, 3); k_third += tmp1 * nodal_area; // eq 107c tmp.mul(t_outer_t, Aj); tmp1.mul(A, tmp); tmp1 *= contravariant_metric_tensor(alpha, theta) * contravariant_metric_tensor(beta, gamma); tmp1 *= mu * p_n * tangential_traction[alpha]; tmp1 /= traction_norm; Matrix tmp2(spatial_dimension, spatial_dimension * nb_nodes_per_contact); Matrix tmp3(nb_nodes_per_contact * spatial_dimension, spatial_dimension * nb_nodes_per_contact); tmp2.mul(t_outer_t, A); tmp3.mul(Aj, tmp2); tmp3 *= contravariant_metric_tensor(alpha, gamma) * contravariant_metric_tensor(beta, theta); tmp3 *= mu * p_n * tangential_traction[alpha]; tmp3 /= traction_norm; k_fourth += (tmp1 + tmp3) * nodal_area; } } } } stiffness += k_third + k_fourth - k_first - k_second; } /* -------------------------------------------------------------------------- */ void ResolutionPenalty::beforeSolveStep() {} /* -------------------------------------------------------------------------- */ void ResolutionPenalty::afterSolveStep(__attribute__((unused)) bool converged) { /*auto method = model.getAnalysisMethod(); if (method == _explicit_lumped_mass) { return ; } auto & K = const_cast(model.getDOFManager().getMatrix("K")); auto k_min = K.min(); auto roundoff_error = 1e-17; const auto blocked_dofs = model.getDOFManager().getBlockedDOFs("displacement"); Real nb_unknowns = 0; for (auto & bld : make_view(blocked_dofs)) { if (not bld) nb_unknowns++; } auto max_epsilon_n = k_min / sqrt(nb_unknowns * roundoff_error); if (epsilon_n > max_epsilon_n) epsilon_n = max_epsilon_n;*/ } INSTANTIATE_RESOLUTION(penalty_linear, ResolutionPenalty); } // namespace akantu diff --git a/src/model/contact_mechanics/resolutions/resolution_penalty_quadratic.cc b/src/model/contact_mechanics/resolutions/resolution_penalty_quadratic.cc index c914ff133..832187d8d 100644 --- a/src/model/contact_mechanics/resolutions/resolution_penalty_quadratic.cc +++ b/src/model/contact_mechanics/resolutions/resolution_penalty_quadratic.cc @@ -1,886 +1,852 @@ /** * @file resolution_penalty_quadratic.cc * * @author Mohit Pundir * * @date creation: Sun Aug 2 2020 * @date last modification: Sun Aug 2 2020 * * @brief Specialization of the resolution class for the quadratic penalty * method * * @section LICENSE * * Copyright (©) 2010-2018 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 "resolution_penalty_quadratic.hh" +#include "element_class_helper.hh" +/* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ ResolutionPenaltyQuadratic::ResolutionPenaltyQuadratic( ContactMechanicsModel & model, const ID & id) : Parent(model, id) { AKANTU_DEBUG_IN(); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ResolutionPenaltyQuadratic::initialize() {} /* -------------------------------------------------------------------------- */ Real ResolutionPenaltyQuadratic::computeNormalTraction(Real & gap) { return epsilon_n * (macaulay(gap) * macaulay(gap) + macaulay(gap)); } /* -------------------------------------------------------------------------- */ void ResolutionPenaltyQuadratic::computeNormalForce( const ContactElement & element, Vector & force) { force.zero(); auto & gaps = model.getGaps(); auto & projections = model.getProjections(); auto & normals = model.getNormals(); auto surface_dimension = spatial_dimension - 1; Real gap(gaps.begin()[element.slave]); Vector normal(normals.begin(spatial_dimension)[element.slave]); Vector projection(projections.begin(surface_dimension)[element.slave]); auto & nodal_area = const_cast &>(model.getNodalArea()); // compute normal traction Real p_n = computeNormalTraction(gap); p_n *= nodal_area[element.slave]; UInt nb_nodes_per_contact = element.getNbNodes(); Matrix shape_matric(spatial_dimension, spatial_dimension * nb_nodes_per_contact); ResolutionUtils::computeShapeFunctionMatric(element, projection, shape_matric); force.mul(shape_matric, normal, p_n); } /* -------------------------------------------------------------------------- */ void ResolutionPenaltyQuadratic::computeTangentialForce( const ContactElement & element, Vector & force) { if (mu == 0) return; force.zero(); UInt surface_dimension = spatial_dimension - 1; // compute tangents auto & projections = model.getProjections(); Vector projection(projections.begin(surface_dimension)[element.slave]); auto & normals = model.getNormals(); Vector normal(normals.begin(spatial_dimension)[element.slave]); auto & tangents = model.getTangents(); Matrix covariant_basis( tangents.begin(surface_dimension, spatial_dimension)[element.slave]); // check for no-contact to contact condition // need a better way to check if new node added is not presnt in the // previous master elemets auto & previous_master_elements = model.getPreviousMasterElements(); if (element.slave >= previous_master_elements.size()) return; auto & previous_element = previous_master_elements[element.slave]; if (previous_element.type == _not_defined) return; // compute tangential traction using return map algorithm auto & tangential_tractions = model.getTangentialTractions(); Vector tangential_traction( tangential_tractions.begin(surface_dimension)[element.slave]); this->computeTangentialTraction(element, covariant_basis, tangential_traction); UInt nb_nodes_per_contact = element.getNbNodes(); Matrix shape_matric(spatial_dimension, spatial_dimension * nb_nodes_per_contact); ResolutionUtils::computeShapeFunctionMatric(element, projection, shape_matric); auto contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(covariant_basis); auto & nodal_area = const_cast &>(model.getNodalArea()); for (auto && values1 : enumerate(covariant_basis.transpose())) { auto & alpha = std::get<0>(values1); auto & tangent_alpha = std::get<1>(values1); for (auto && values2 : enumerate(tangential_traction)) { auto & beta = std::get<0>(values2); auto & traction_beta = std::get<1>(values2); Vector tmp(force.size()); tmp.mul(shape_matric, tangent_alpha, traction_beta); tmp *= contravariant_metric_tensor(alpha, beta) * nodal_area[element.slave]; force += tmp; } } // compute first variation of natural coordinate /*auto & gaps = model.getGaps(); auto & gap = gaps.begin()[element.slave]; auto nb_nodes = element.getNbNodes(); Array delta_xi(nb_nodes * spatial_dimension, surface_dimension); ResolutionUtils::firstVariationNaturalCoordinate(element, covariant_basis, projection, normal, gap, delta_xi); // compute tangential force auto & nodal_area = const_cast &>(model.getNodalArea()); for (auto && values : zip(tangential_traction, make_view(delta_xi, delta_xi.size()))) { auto & traction_alpha = std::get<0>(values); auto & delta_xi_alpha = std::get<1>(values); force += delta_xi_alpha * traction_alpha * nodal_area[element.slave]; }*/ } /* -------------------------------------------------------------------------- */ void ResolutionPenaltyQuadratic::computeTangentialTraction( const ContactElement & element, const Matrix & covariant_basis, Vector & traction_tangential) { UInt surface_dimension = spatial_dimension - 1; auto & gaps = model.getGaps(); auto & gap = gaps.begin()[element.slave]; // Return map algorithm is employed // compute trial traction Vector traction_trial(surface_dimension); this->computeTrialTangentialTraction(element, covariant_basis, traction_trial); // compute norm of trial traction Real traction_trial_norm = 0; auto contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(covariant_basis); for (auto i : arange(surface_dimension)) { for (auto j : arange(surface_dimension)) { traction_trial_norm += traction_trial[i] * traction_trial[j] * contravariant_metric_tensor(i, j); } } traction_trial_norm = sqrt(traction_trial_norm); // check stick or slip condition auto & contact_state = model.getContactState(); auto & state = contact_state.begin()[element.slave]; Real p_n = computeNormalTraction(gap); bool stick = (traction_trial_norm <= mu * p_n) ? true : false; if (stick) { state = ContactState::_stick; computeStickTangentialTraction(element, traction_trial, traction_tangential); } else { state = ContactState::_slip; computeSlipTangentialTraction(element, covariant_basis, traction_trial, traction_tangential); } } /* -------------------------------------------------------------------------- */ void ResolutionPenaltyQuadratic::computeTrialTangentialTraction( const ContactElement & element, const Matrix & covariant_basis, Vector & traction) { UInt surface_dimension = spatial_dimension - 1; auto & projections = model.getProjections(); Vector current_projection( projections.begin(surface_dimension)[element.slave]); auto & previous_projections = model.getPreviousProjections(); Vector previous_projection( previous_projections.begin(surface_dimension)[element.slave]); // method from Laursen et. al. /*auto covariant_metric_tensor = GeometryUtils::covariantMetricTensor(covariant_basis); auto increment_projection = current_projection - previous_projection; traction.mul(covariant_metric_tensor, increment_projection, epsilon_t); auto & previous_tangential_tractions = model.getPreviousTangentialTractions(); Vector previous_traction(previous_tangential_tractions.begin(surface_dimension)[element.slave]); traction = previous_traction + traction;*/ // method from Schweizerhof auto covariant_metric_tensor = GeometryUtils::covariantMetricTensor(covariant_basis); auto & previous_tangential_tractions = model.getPreviousTangentialTractions(); Vector previous_traction( previous_tangential_tractions.begin(surface_dimension)[element.slave]); auto & previous_tangents = model.getPreviousTangents(); Matrix previous_covariant_basis(previous_tangents.begin( surface_dimension, spatial_dimension)[element.slave]); auto previous_contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(previous_covariant_basis); auto current_tangent = covariant_basis.transpose(); auto previous_tangent = previous_covariant_basis.transpose(); for (auto alpha : arange(surface_dimension)) { Vector tangent_alpha(current_tangent(alpha)); for (auto gamma : arange(surface_dimension)) { for (auto beta : arange(surface_dimension)) { Vector tangent_beta(previous_tangent(beta)); auto t_alpha_t_beta = tangent_beta.dot(tangent_alpha); traction[alpha] += previous_traction[gamma] * previous_contravariant_metric_tensor(gamma, beta) * t_alpha_t_beta; } } } auto & previous_master_elements = model.getPreviousMasterElements(); auto & previous_element = previous_master_elements[element.slave]; Vector previous_real_projection(spatial_dimension); GeometryUtils::realProjection( model.getMesh(), model.getContactDetector().getPositions(), previous_element, previous_projection, previous_real_projection); Vector current_real_projection(spatial_dimension); GeometryUtils::realProjection( model.getMesh(), model.getContactDetector().getPositions(), element.master, current_projection, current_real_projection); auto increment_real = current_real_projection - previous_real_projection; Vector increment_xi(surface_dimension); auto contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(covariant_basis); // increment in natural coordinate for (auto beta : arange(surface_dimension)) { for (auto gamma : arange(surface_dimension)) { auto temp = increment_real.dot(current_tangent(gamma)); temp *= contravariant_metric_tensor(beta, gamma); increment_xi[beta] += temp; } } Vector temp(surface_dimension); temp.mul(covariant_metric_tensor, increment_xi, epsilon_t); traction -= temp; } /* -------------------------------------------------------------------------- */ void ResolutionPenaltyQuadratic::computeStickTangentialTraction( const ContactElement & /*element*/, Vector & traction_trial, Vector & traction_tangential) { traction_tangential = traction_trial; } /* -------------------------------------------------------------------------- */ void ResolutionPenaltyQuadratic::computeSlipTangentialTraction( const ContactElement & element, const Matrix & covariant_basis, Vector & traction_trial, Vector & traction_tangential) { UInt surface_dimension = spatial_dimension - 1; auto & gaps = model.getGaps(); auto & gap = gaps.begin()[element.slave]; // compute norm of trial traction Real traction_trial_norm = 0; auto contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(covariant_basis); for (auto i : arange(surface_dimension)) { for (auto j : arange(surface_dimension)) { traction_trial_norm += traction_trial[i] * traction_trial[j] * contravariant_metric_tensor(i, j); } } traction_trial_norm = sqrt(traction_trial_norm); auto slip_direction = traction_trial; slip_direction /= traction_trial_norm; Real p_n = computeNormalTraction(gap); traction_tangential = slip_direction; traction_tangential *= mu * p_n; } /* -------------------------------------------------------------------------- */ void ResolutionPenaltyQuadratic::computeNormalModuli( const ContactElement & element, Matrix & stiffness) { auto surface_dimension = spatial_dimension - 1; auto & gaps = model.getGaps(); Real gap(gaps.begin()[element.slave]); auto & projections = model.getProjections(); Vector projection(projections.begin(surface_dimension)[element.slave]); auto & nodal_areas = model.getNodalArea(); auto & nodal_area = nodal_areas.begin()[element.slave]; auto & normals = model.getNormals(); Vector normal(normals.begin(spatial_dimension)[element.slave]); - auto & mesh = model.getMesh(); - // method from Schweizerhof and A. Konyukhov, K. Schweizerhof // DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3 // construct A matrix const ElementType & type = element.master.type; - UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); - Vector shapes(nb_nodes_per_element); -#define GET_SHAPE_NATURAL(type) \ - ElementClass::computeShapes(projection, shapes) - AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_NATURAL); -#undef GET_SHAPE_NATURAL + auto && shapes = ElementClassHelper<_ek_regular>::getN(projection, type); UInt nb_nodes_per_contact = element.getNbNodes(); Matrix A(spatial_dimension, spatial_dimension * nb_nodes_per_contact); for (auto i : arange(nb_nodes_per_contact)) { for (auto j : arange(spatial_dimension)) { if (i == 0) { A(j, i * spatial_dimension + j) = 1; continue; } A(j, i * spatial_dimension + j) = -shapes[i - 1]; } } // construct the main part of normal matrix Matrix k_main(nb_nodes_per_contact * spatial_dimension, nb_nodes_per_contact * spatial_dimension); Matrix n_outer_n(spatial_dimension, spatial_dimension); Matrix mat_n(normal.storage(), normal.size(), 1.); n_outer_n.mul(mat_n, mat_n); Matrix tmp(spatial_dimension, spatial_dimension * nb_nodes_per_contact); tmp.mul(n_outer_n, A); k_main.mul(A, tmp); k_main *= epsilon_n * heaviside(gap) * (2 * gap + 1) * nodal_area; // construct the rotational part of the normal matrix auto & tangents = model.getTangents(); Matrix covariant_basis( tangents.begin(surface_dimension, spatial_dimension)[element.slave]); auto contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(covariant_basis); // computing shape derivatives - Matrix shape_derivatives(surface_dimension, - Mesh::getNbNodesPerElement(type)); - -#define GET_SHAPE_DERIVATIVES_NATURAL(type) \ - ElementClass::computeDNDS(projection, shape_derivatives) - AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); -#undef GET_SHAPE_DERIVATIVES_NATURAL + auto && shape_derivatives = + ElementClassHelper<_ek_regular>::getDNDS(projection, type); // consists of 2 rotational parts Matrix k_rot1(nb_nodes_per_contact * spatial_dimension, nb_nodes_per_contact * spatial_dimension); Matrix k_rot2(nb_nodes_per_contact * spatial_dimension, nb_nodes_per_contact * spatial_dimension); Matrix Aj(spatial_dimension, spatial_dimension * nb_nodes_per_contact); auto construct_Aj = [&](auto && dnds) { for (auto i : arange(nb_nodes_per_contact)) { for (auto j : arange(spatial_dimension)) { if (i == 0) { Aj(j, i * spatial_dimension + j) = 0; continue; } Aj(j, i * spatial_dimension + j) = dnds(i - 1); } } }; for (auto && values1 : enumerate(covariant_basis.transpose())) { auto & alpha = std::get<0>(values1); auto & tangent = std::get<1>(values1); Matrix n_outer_t(spatial_dimension, spatial_dimension); Matrix mat_t(tangent.storage(), tangent.size(), 1.); n_outer_t.mul(mat_n, mat_t); Matrix t_outer_n(spatial_dimension, spatial_dimension); t_outer_n.mul(mat_t, mat_n); for (auto && values2 : enumerate(shape_derivatives.transpose())) { auto & beta = std::get<0>(values2); auto & dnds = std::get<1>(values2); // construct Aj from shape function wrt to jth natural // coordinate construct_Aj(dnds); Matrix tmp(spatial_dimension, spatial_dimension * nb_nodes_per_contact); Matrix tmp1(nb_nodes_per_contact * spatial_dimension, spatial_dimension * nb_nodes_per_contact); tmp.mul(n_outer_t, A); tmp1.mul(Aj, tmp); tmp1 *= contravariant_metric_tensor(alpha, beta); k_rot1 += tmp1; tmp.mul(t_outer_n, Aj); tmp1.mul(A, tmp); tmp1 *= contravariant_metric_tensor(alpha, beta); k_rot2 += tmp1; } } k_rot1 *= -epsilon_n * heaviside(gap) * (gap * gap + gap) * nodal_area; k_rot2 *= -epsilon_n * heaviside(gap) * (gap * gap + gap) * nodal_area; stiffness += k_main + k_rot1 + k_rot2; } /* -------------------------------------------------------------------------- */ void ResolutionPenaltyQuadratic::computeTangentialModuli( const ContactElement & element, Matrix & stiffness) { if (mu == 0) return; stiffness.zero(); auto & contact_state = model.getContactState(); UInt state = contact_state.begin()[element.slave]; switch (state) { case ContactState::_stick: { computeStickModuli(element, stiffness); break; } case ContactState::_slip: { computeSlipModuli(element, stiffness); break; } default: break; } } /* -------------------------------------------------------------------------- */ void ResolutionPenaltyQuadratic::computeStickModuli( const ContactElement & element, Matrix & stiffness) { auto surface_dimension = spatial_dimension - 1; auto & projections = model.getProjections(); Vector projection(projections.begin(surface_dimension)[element.slave]); auto & nodal_areas = model.getNodalArea(); auto & nodal_area = nodal_areas.begin()[element.slave]; - auto & mesh = model.getMesh(); - // method from Schweizerhof and A. Konyukhov, K. Schweizerhof // DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3 // construct A matrix const ElementType & type = element.master.type; - UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); - Vector shapes(nb_nodes_per_element); - -#define GET_SHAPE_NATURAL(type) \ - ElementClass::computeShapes(projection, shapes) - AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_NATURAL); -#undef GET_SHAPE_NATURAL + auto && shapes = ElementClassHelper<_ek_regular>::getN(projection, type); UInt nb_nodes_per_contact = element.getNbNodes(); Matrix A(spatial_dimension, spatial_dimension * nb_nodes_per_contact); for (auto i : arange(nb_nodes_per_contact)) { for (auto j : arange(spatial_dimension)) { if (i == 0) { A(j, i * spatial_dimension + j) = 1; continue; } A(j, i * spatial_dimension + j) = -shapes[i - 1]; } } // computing shape derivatives - Matrix shape_derivatives(surface_dimension, - Mesh::getNbNodesPerElement(type)); - -#define GET_SHAPE_DERIVATIVES_NATURAL(type) \ - ElementClass::computeDNDS(projection, shape_derivatives) - AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); -#undef GET_SHAPE_DERIVATIVES_NATURAL + auto && shape_derivatives = + ElementClassHelper<_ek_regular>::getDNDS(projection, type); Matrix Aj(spatial_dimension, spatial_dimension * nb_nodes_per_contact); auto construct_Aj = [&](auto && dnds) { for (auto i : arange(nb_nodes_per_contact)) { for (auto j : arange(spatial_dimension)) { if (i == 0) { Aj(j, i * spatial_dimension + j) = 0; continue; } Aj(j, i * spatial_dimension + j) = dnds(i - 1); } } }; // tangents should have been calculated in normal modulii auto & tangents = model.getTangents(); Matrix covariant_basis( tangents.begin(surface_dimension, spatial_dimension)[element.slave]); auto contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(covariant_basis); // construct 1st part of the stick modulii Matrix k_main(nb_nodes_per_contact * spatial_dimension, nb_nodes_per_contact * spatial_dimension); for (auto && values1 : enumerate(covariant_basis.transpose())) { auto & alpha = std::get<0>(values1); auto & tangent_alpha = std::get<1>(values1); Matrix t_outer_t(spatial_dimension, spatial_dimension); Matrix mat_t_alpha(tangent_alpha.storage(), tangent_alpha.size(), 1.); for (auto && values2 : enumerate(covariant_basis.transpose())) { auto & beta = std::get<0>(values2); auto & tangent_beta = std::get<1>(values2); Matrix mat_t_beta(tangent_beta.storage(), tangent_beta.size(), 1.); t_outer_t.mul(mat_t_alpha, mat_t_beta); Matrix tmp(spatial_dimension, spatial_dimension * nb_nodes_per_contact); Matrix tmp1(nb_nodes_per_contact * spatial_dimension, spatial_dimension * nb_nodes_per_contact); tmp.mul(t_outer_t, A); tmp1.mul(A, tmp); tmp1 *= contravariant_metric_tensor(alpha, beta); k_main += tmp1; } } k_main *= -epsilon_t; // construct 2nd part of the stick modulii auto & tangential_tractions = model.getTangentialTractions(); Vector tangential_traction( tangential_tractions.begin(surface_dimension)[element.slave]); Matrix k_second(nb_nodes_per_contact * spatial_dimension, nb_nodes_per_contact * spatial_dimension); for (auto alpha : arange(surface_dimension)) { Matrix k_sum(nb_nodes_per_contact * spatial_dimension, nb_nodes_per_contact * spatial_dimension); for (auto && values1 : enumerate(shape_derivatives.transpose())) { auto & beta = std::get<0>(values1); auto & dnds = std::get<1>(values1); // construct Aj from shape function wrt to jth natural // coordinate construct_Aj(dnds); for (auto && values2 : enumerate(covariant_basis.transpose())) { auto & gamma = std::get<0>(values2); auto & tangent_gamma = std::get<1>(values2); Matrix t_outer_t(spatial_dimension, spatial_dimension); Matrix mat_t_gamma(tangent_gamma.storage(), tangent_gamma.size(), 1.); for (auto && values3 : enumerate(covariant_basis.transpose())) { auto & theta = std::get<0>(values3); auto & tangent_theta = std::get<1>(values3); Matrix mat_t_theta(tangent_theta.storage(), tangent_theta.size(), 1.); t_outer_t.mul(mat_t_gamma, mat_t_theta); Matrix tmp(spatial_dimension, spatial_dimension * nb_nodes_per_contact); Matrix tmp1(nb_nodes_per_contact * spatial_dimension, spatial_dimension * nb_nodes_per_contact); tmp.mul(t_outer_t, Aj); tmp1.mul(A, tmp); tmp1 *= contravariant_metric_tensor(alpha, theta) * contravariant_metric_tensor(beta, gamma); Matrix tmp2(spatial_dimension, spatial_dimension * nb_nodes_per_contact); Matrix tmp3(nb_nodes_per_contact * spatial_dimension, spatial_dimension * nb_nodes_per_contact); tmp2.mul(t_outer_t, A); tmp3.mul(Aj, tmp2); tmp3 *= contravariant_metric_tensor(alpha, gamma) * contravariant_metric_tensor(beta, theta); k_sum += tmp1 + tmp3; } } } k_second += tangential_traction[alpha] * k_sum; } stiffness += k_main * nodal_area - k_second * nodal_area; } /* -------------------------------------------------------------------------- */ void ResolutionPenaltyQuadratic::computeSlipModuli( const ContactElement & element, Matrix & stiffness) { auto surface_dimension = spatial_dimension - 1; auto & gaps = model.getGaps(); Real gap(gaps.begin()[element.slave]); auto & nodal_areas = model.getNodalArea(); auto & nodal_area = nodal_areas.begin()[element.slave]; // compute normal traction Real p_n = computeNormalTraction(gap); auto & projections = model.getProjections(); Vector projection(projections.begin(surface_dimension)[element.slave]); auto & normals = model.getNormals(); Vector normal(normals.begin(spatial_dimension)[element.slave]); // restructure normal as a matrix for an outer product Matrix mat_n(normal.storage(), normal.size(), 1.); auto & mesh = model.getMesh(); // method from Schweizerhof and A. Konyukhov, K. Schweizerhof // DOI 10.1007/s00466-004-0616-7 and DOI 10.1007/s00466-003-0515-3 // construct A matrix const ElementType & type = element.master.type; - UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); - Vector shapes(nb_nodes_per_element); - -#define GET_SHAPE_NATURAL(type) \ - ElementClass::computeShapes(projection, shapes) - AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_NATURAL); -#undef GET_SHAPE_NATURAL + auto && shapes = ElementClassHelper<_ek_regular>::getN(projection, type); UInt nb_nodes_per_contact = element.getNbNodes(); Matrix A(spatial_dimension, spatial_dimension * nb_nodes_per_contact); for (auto i : arange(nb_nodes_per_contact)) { for (auto j : arange(spatial_dimension)) { if (i == 0) { A(j, i * spatial_dimension + j) = 1; continue; } A(j, i * spatial_dimension + j) = -shapes[i - 1]; } } // computing shape derivatives - Matrix shape_derivatives(surface_dimension, - Mesh::getNbNodesPerElement(type)); - -#define GET_SHAPE_DERIVATIVES_NATURAL(type) \ - ElementClass::computeDNDS(projection, shape_derivatives) - AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_NATURAL); -#undef GET_SHAPE_DERIVATIVES_NATURAL + auto && shape_derivatives = + ElementClassHelper<_ek_regular>::getDNDS(projection, type); Matrix Aj(spatial_dimension, spatial_dimension * nb_nodes_per_contact); auto construct_Aj = [&](auto && dnds) { for (auto i : arange(nb_nodes_per_contact)) { for (auto j : arange(spatial_dimension)) { if (i == 0) { Aj(j, i * spatial_dimension + j) = 0; continue; } Aj(j, i * spatial_dimension + j) = dnds(i - 1); } } }; // tangents should have been calculated in normal modulii auto & tangents = model.getTangents(); Matrix covariant_basis( tangents.begin(surface_dimension, spatial_dimension)[element.slave]); auto & tangential_tractions = model.getTangentialTractions(); Vector tangential_traction( tangential_tractions.begin(surface_dimension)[element.slave]); // compute norm of trial traction Real traction_norm = 0; auto contravariant_metric_tensor = GeometryUtils::contravariantMetricTensor(covariant_basis); for (auto i : arange(surface_dimension)) { for (auto j : arange(surface_dimension)) { traction_norm += tangential_traction[i] * tangential_traction[j] * contravariant_metric_tensor(i, j); } } traction_norm = sqrt(traction_norm); // construct four parts of stick modulii (eq 107,107a-c) Matrix k_first(nb_nodes_per_contact * spatial_dimension, nb_nodes_per_contact * spatial_dimension); Matrix k_second(nb_nodes_per_contact * spatial_dimension, nb_nodes_per_contact * spatial_dimension); Matrix k_third(nb_nodes_per_contact * spatial_dimension, nb_nodes_per_contact * spatial_dimension); Matrix k_fourth(nb_nodes_per_contact * spatial_dimension, nb_nodes_per_contact * spatial_dimension); for (auto && values1 : enumerate(covariant_basis.transpose())) { auto & alpha = std::get<0>(values1); auto & tangent_alpha = std::get<1>(values1); Matrix mat_t_alpha(tangent_alpha.storage(), tangent_alpha.size(), 1.); Matrix t_outer_n(spatial_dimension, spatial_dimension); Matrix t_outer_t(spatial_dimension, spatial_dimension); for (auto && values2 : zip(arange(surface_dimension), covariant_basis.transpose(), shape_derivatives.transpose())) { auto & beta = std::get<0>(values2); auto & tangent_beta = std::get<1>(values2); auto & dnds = std::get<2>(values2); // construct Aj from shape function wrt to jth natural // coordinate construct_Aj(dnds); // eq 107 Matrix mat_t_beta(tangent_beta.storage(), tangent_beta.size(), 1.); t_outer_n.mul(mat_t_beta, mat_n); Matrix tmp(spatial_dimension, spatial_dimension * nb_nodes_per_contact); Matrix tmp1(nb_nodes_per_contact * spatial_dimension, spatial_dimension * nb_nodes_per_contact); tmp.mul(t_outer_n, A); tmp1.mul(A, tmp); tmp1 *= epsilon_n * mu * tangential_traction[alpha] * contravariant_metric_tensor(alpha, beta); tmp1 /= traction_norm; k_first += tmp1 * nodal_area; // eq 107a t_outer_t.mul(mat_t_alpha, mat_t_beta); tmp.mul(t_outer_t, A); tmp1.mul(A, tmp); tmp1 *= epsilon_t * mu * p_n * contravariant_metric_tensor(alpha, beta); tmp1 /= traction_norm; k_second += tmp1 * nodal_area; for (auto && values3 : enumerate(covariant_basis.transpose())) { auto & gamma = std::get<0>(values3); auto & tangent_gamma = std::get<1>(values3); Matrix mat_t_gamma(tangent_gamma.storage(), tangent_gamma.size(), 1.); for (auto && values4 : enumerate(covariant_basis.transpose())) { auto & theta = std::get<0>(values4); auto & tangent_theta = std::get<1>(values4); Matrix mat_t_theta(tangent_theta.storage(), tangent_theta.size(), 1.); t_outer_t.mul(mat_t_gamma, mat_t_theta); // eq 107b tmp.mul(t_outer_t, A); tmp1.mul(A, tmp); tmp1 *= epsilon_t * mu * p_n * tangential_traction[alpha] * tangential_traction[beta]; tmp1 *= contravariant_metric_tensor(alpha, gamma) * contravariant_metric_tensor(beta, theta); tmp1 /= pow(traction_norm, 3); k_third += tmp1 * nodal_area; // eq 107c tmp.mul(t_outer_t, Aj); tmp1.mul(A, tmp); tmp1 *= contravariant_metric_tensor(alpha, theta) * contravariant_metric_tensor(beta, gamma); tmp1 *= mu * p_n * tangential_traction[alpha]; tmp1 /= traction_norm; Matrix tmp2(spatial_dimension, spatial_dimension * nb_nodes_per_contact); Matrix tmp3(nb_nodes_per_contact * spatial_dimension, spatial_dimension * nb_nodes_per_contact); tmp2.mul(t_outer_t, A); tmp3.mul(Aj, tmp2); tmp3 *= contravariant_metric_tensor(alpha, gamma) * contravariant_metric_tensor(beta, theta); tmp3 *= mu * p_n * tangential_traction[alpha]; tmp3 /= traction_norm; k_fourth += (tmp1 + tmp3) * nodal_area; } } } } stiffness += k_third + k_fourth - k_first - k_second; } /* -------------------------------------------------------------------------- */ void ResolutionPenaltyQuadratic::beforeSolveStep() {} /* -------------------------------------------------------------------------- */ void ResolutionPenaltyQuadratic::afterSolveStep( __attribute__((unused)) bool converged) {} INSTANTIATE_RESOLUTION(penalty_quadratic, ResolutionPenaltyQuadratic); } // namespace akantu