diff --git a/packages/05_igfem.cmake b/packages/05_igfem.cmake index e4a3ab330..e0831bdb2 100644 --- a/packages/05_igfem.cmake +++ b/packages/05_igfem.cmake @@ -1,37 +1,43 @@ #=============================================================================== # @file igfem.cmake # # @author Aurelia Cuba Ramos # @author Nicolas Richart # # @date Fri Mai 23 18:19:15 2011 # # @brief package description for interface-enriched generalized IGFEM # # @section LICENSE # # Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # #=============================================================================== option(AKANTU_IGFEM "Use Interface-enriched generalized FEM" OFF) set(AKANTU_IGFEM_FILES fem/element_class_igfem.hh fem/shape_igfem.hh fem/shape_igfem_inline_impl.cc + fem/igfem_element.hh ) + +set(AKANTU_IGFEM_TESTS + test_igfem_integrate +) + diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh index 2eb305b0c..cfef514f9 100644 --- a/src/common/aka_common.hh +++ b/src/common/aka_common.hh @@ -1,583 +1,591 @@ /** * @file aka_common.hh * * @author Nicolas Richart * * @date Mon Jun 14 19:12:20 2010 * * @brief common type descriptions for akantu * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * * @section DESCRIPTION * * All common things to be included in the projects files * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COMMON_HH__ #define __AKANTU_COMMON_HH__ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #define __BEGIN_AKANTU__ namespace akantu { #define __END_AKANTU__ }; /* -------------------------------------------------------------------------- */ #include "aka_config.hh" #include "aka_error.hh" #include "aka_safe_enum.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* Common types */ /* -------------------------------------------------------------------------- */ typedef double Real; typedef unsigned int UInt; typedef unsigned long long UInt64; typedef signed int Int; typedef std::string ID; static const Real UINT_INIT_VALUE = 0; #ifdef AKANTU_NDEBUG static const Real REAL_INIT_VALUE = 0; #else static const Real REAL_INIT_VALUE = std::numeric_limits::quiet_NaN(); #endif /* -------------------------------------------------------------------------- */ /* Memory types */ /* -------------------------------------------------------------------------- */ typedef UInt MemoryID; /* -------------------------------------------------------------------------- */ /* Mesh/FEM/Model types */ /* -------------------------------------------------------------------------- */ typedef std::string Surface; typedef std::pair SurfacePair; typedef std::list< SurfacePair > SurfacePairList; /* -------------------------------------------------------------------------- */ extern const UInt _all_dimensions; /// @boost sequence of element to loop on in global tasks #define AKANTU_REGULAR_ELEMENT_TYPE \ (_point_1) \ (_segment_2) \ (_segment_3) \ (_triangle_3) \ (_triangle_6) \ (_quadrangle_4) \ (_quadrangle_8) \ (_tetrahedron_4) \ (_tetrahedron_10) \ (_pentahedron_6) \ (_hexahedron_8) #if defined(AKANTU_STRUCTURAL_MECHANICS) #define AKANTU_STRUCTURAL_ELEMENT_TYPE \ (_bernoulli_beam_2) \ (_bernoulli_beam_3) #else #define AKANTU_STRUCTURAL_ELEMENT_TYPE #endif #if defined(AKANTU_COHESIVE_ELEMENT) # define AKANTU_COHESIVE_ELEMENT_TYPE \ (_cohesive_2d_4) \ (_cohesive_2d_6) \ (_cohesive_1d_2) \ (_cohesive_3d_6) \ (_cohesive_3d_12) #else # define AKANTU_COHESIVE_ELEMENT_TYPE #endif #if defined(AKANTU_IGFEM) #define AKANTU_IGFEM_ELEMENT_TYPE \ (_igfem_segment_2) \ (_igfem_triangle_3) #else #define AKANTU_IGFEM_ELEMENT_TYPE #endif #define AKANTU_ALL_ELEMENT_TYPE \ AKANTU_REGULAR_ELEMENT_TYPE \ AKANTU_COHESIVE_ELEMENT_TYPE \ AKANTU_STRUCTURAL_ELEMENT_TYPE \ AKANTU_IGFEM_ELEMENT_TYPE #define AKANTU_NOT_STRUCTURAL_ELEMENT_TYPE \ AKANTU_REGULAR_ELEMENT_TYPE \ - AKANTU_COHESIVE_ELEMENT_TYPE + AKANTU_COHESIVE_ELEMENT_TYPE \ + AKANTU_IGFEM_ELEMENT_TYPE /// @enum ElementType type of elements enum ElementType { _not_defined, _point_1, _segment_2, ///< first order segment _segment_3, ///< second order segment _triangle_3, ///< first order triangle _triangle_6, ///< second order triangle _tetrahedron_4, ///< first order tetrahedron _tetrahedron_10, ///< second order tetrahedron _quadrangle_4, ///< first order quadrangle _quadrangle_8, ///< second order quadrangle _hexahedron_8, ///< first order hexahedron _pentahedron_6, ///< first order pentahedron _bernoulli_beam_2, ///< Bernoulli beam 2D _bernoulli_beam_3, ///< Bernoulli beam 3D #if defined(AKANTU_COHESIVE_ELEMENT) _cohesive_2d_4, ///< first order 2D cohesive _cohesive_2d_6, ///< second order 2D cohesive _cohesive_1d_2, ///< first order 1D cohesive _cohesive_3d_6, ///< first order 3D cohesive _cohesive_3d_12, ///< second order 3D cohesive #endif -#if defined(AKANTU_IGFEM_ELEMENT) +#if defined(AKANTU_IGFEM) _igfem_segment_2, ///< first order segment _igfem_triangle_3, ///< first order triangle #endif _max_element_type }; /// @enum GeometricalType type of element potentially contained in a Mesh enum GeometricalType { _gt_point, ///< point @remark only for some algorithm to be generic like mesh partitioning */ _gt_segment_2, ///< 2 nodes segment _gt_segment_3, ///< 3 nodes segment _gt_triangle_3, ///< 3 nodes triangle _gt_triangle_6, ///< 6 nodes triangle _gt_quadrangle_4, ///< 4 nodes quadrangle _gt_quadrangle_8, ///< 8 nodes quadrangle _gt_tetrahedron_4, ///< 4 nodes tetrahedron _gt_tetrahedron_10, ///< 10 nodes tetrahedron _gt_hexahedron_8, ///< 8 nodes hexahedron _gt_pentahedron_6, ///< 6 nodes pentahedron #if defined(AKANTU_COHESIVE_ELEMENT) _gt_cohesive_2d_4, ///< 4 nodes 2D cohesive _gt_cohesive_2d_6, ///< 6 nodes 2D cohesive _gt_cohesive_1d_2, ///< 2 nodes 1D cohesive _gt_cohesive_3d_6, ///< 6 nodes 3D cohesive _gt_cohesive_3d_12, ///< 12 nodes 3D cohesive #endif _gt_not_defined }; /// @enum InterpolationType type of elements enum InterpolationType { _itp_lagrange_point_1, ///< zeroth (!) order lagrangian point (for compatibility purposes) _itp_lagrange_segment_2, ///< first order lagrangian segment _itp_lagrange_segment_3, ///< second order lagrangian segment _itp_lagrange_triangle_3, ///< first order lagrangian triangle _itp_lagrange_triangle_6, ///< second order lagrangian triangle _itp_lagrange_quadrangle_4, ///< first order lagrangian quadrangle _itp_serendip_quadrangle_8, /**< second order serendipian quadrangle @remark used insted of the 9 node lagrangian element */ _itp_lagrange_tetrahedron_4, ///< first order lagrangian tetrahedron _itp_lagrange_tetrahedron_10, ///< second order lagrangian tetrahedron _itp_lagrange_hexahedron_8, ///< first order lagrangian hexahedron _itp_lagrange_pentahedron_6, ///< first order lagrangian pentahedron _itp_bernoulli_beam, ///< Bernoulli beam _itp_not_defined }; //! standard output stream operator for ElementType inline std::ostream & operator <<(std::ostream & stream, ElementType type); enum ElementKind { _ek_regular, _ek_cohesive, _ek_structural, _ek_igfem, _ek_not_defined }; /// enum MeshIOType type of mesh reader/writer enum MeshIOType { _miot_auto, _miot_gmsh, _miot_diana, _miot_abaqus, }; /// enum AnalysisMethod type of solving method used to solve the equation of motion enum AnalysisMethod { _static, _implicit_dynamic, _explicit_lumped_mass, _explicit_lumped_capacity, _explicit_consistent_mass }; /// enum SolveConvergenceMethod different resolution algorithms enum SolveConvergenceMethod { _scm_newton_raphson_tangent, ///< Newton-Raphson with tangent matrix _scm_newton_raphson_tangent_modified, ///< Newton-Raphson with constant tangent matrix _scm_newton_raphson_tangent_not_computed ///< Newton-Raphson with constant tangent matrix (user has to assemble it) }; /// enum SolveConvergenceCriteria different convergence criteria enum SolveConvergenceCriteria { _scc_residual, ///< Use residual to test the convergence _scc_increment ///< Use increment to test the convergence }; /// enum CohesiveMethod type of insertion of cohesive elements enum CohesiveMethod { _intrinsic, _extrinsic }; /// myfunction(double * position, double * stress/force, double * normal, unsigned int material_id) typedef void (*BoundaryFunction)(double *,double *, double*, unsigned int); /// @enum BoundaryFunctionType type of function passed for boundary conditions enum BoundaryFunctionType { _bft_stress, _bft_traction, _bft_traction_local }; /// @enum SparseMatrixType type of sparse matrix used enum SparseMatrixType { _unsymmetric, _symmetric }; /* -------------------------------------------------------------------------- */ /* Contact */ /* -------------------------------------------------------------------------- */ typedef ID ContactID; typedef ID ContactSearchID; typedef ID ContactNeighborStructureID; enum ContactType { _ct_not_defined = 0, _ct_2d_expli = 1, _ct_3d_expli = 2, _ct_rigid = 3 }; enum ContactSearchType { _cst_not_defined = 0, _cst_2d_expli = 1, _cst_expli = 2 }; enum ContactNeighborStructureType { _cnst_not_defined = 0, _cnst_regular_grid = 1, _cnst_2d_grid = 2 }; /* -------------------------------------------------------------------------- */ /* Friction */ /* -------------------------------------------------------------------------- */ typedef ID FrictionID; /* -------------------------------------------------------------------------- */ /* Ghosts handling */ /* -------------------------------------------------------------------------- */ typedef ID SynchronizerID; /// @enum CommunicatorType type of communication method to use enum CommunicatorType { _communicator_mpi, _communicator_dummy }; /// @enum SynchronizationTag type of synchronizations enum SynchronizationTag { //--- SolidMechanicsModel tags --- _gst_smm_mass, //< synchronization of the SolidMechanicsModel.mass _gst_smm_for_strain, //< synchronization of the SolidMechanicsModel.current_position _gst_smm_boundary, //< synchronization of the boundary, forces, velocities and displacement _gst_smm_uv, //< synchronization of the nodal velocities and displacement _gst_smm_res, //< synchronization of the nodal residual _gst_smm_init_mat, //< synchronization of the data to initialize materials _gst_smm_stress, //< synchronization of the stresses to compute the internal forces _gst_smmc_facets, //< synchronization of facet data to setup facet synch _gst_smmc_facets_conn, //< synchronization of facet global connectivity _gst_smmc_facets_stress, //< synchronization of facets' stress to setup facet synch _gst_smmc_damage, //< synchronization of damage //--- CohesiveElementInserter tags --- _gst_ce_inserter, //< synchronization of global nodes id of newly inserted cohesive elements //--- GroupManager tags --- _gst_gm_clusters, //< synchronization of clusters //--- HeatTransfer tags --- _gst_htm_capacity, //< synchronization of the nodal heat capacity _gst_htm_temperature, //< synchronization of the nodal temperature _gst_htm_gradient_temperature, //< synchronization of the element gradient temperature //--- LevelSet tags --- /// synchronization of the nodal level set value phi _gst_htm_phi, /// synchronization of the element gradient phi _gst_htm_gradient_phi, //--- Material non local --- _gst_mnl_for_average, //< synchronization of data to average in non local material _gst_mnl_weight, //< synchronization of data for the weight computations //--- General tags --- _gst_test, //< Test tag _gst_material_id, //< synchronization of the material ids //--- Contact & Friction --- _gst_cf_nodal, //< synchronization of disp, velo, and current position _gst_cf_incr //< synchronization of increment }; /// standard output stream operator for SynchronizationTag inline std::ostream & operator <<(std::ostream & stream, SynchronizationTag type); /// @enum GhostType type of ghost enum GhostType { _not_ghost, _ghost, _casper // not used but a real cute ghost }; /* -------------------------------------------------------------------------- */ struct GhostType_def { typedef GhostType type; static const type _begin_ = _not_ghost; static const type _end_ = _casper; }; typedef safe_enum ghost_type_t; /// standard output stream operator for GhostType inline std::ostream & operator <<(std::ostream & stream, GhostType type); /// @enum SynchronizerOperation reduce operation that the synchronizer can perform enum SynchronizerOperation { _so_sum, _so_min, _so_max, _so_null }; /* -------------------------------------------------------------------------- */ /* Global defines */ /* -------------------------------------------------------------------------- */ #define AKANTU_MIN_ALLOCATION 2000 #define AKANTU_INDENT " " #define AKANTU_INCLUDE_INLINE_IMPL /* -------------------------------------------------------------------------- */ // int 2 type construct template struct Int2Type { enum { result = d }; }; // type 2 type construct template class Type2Type { typedef T OriginalType; }; /* -------------------------------------------------------------------------- */ template struct is_scalar { enum{ value = false }; }; #define AKANTU_SPECIFY_IS_SCALAR(type) \ template<> \ struct is_scalar { \ enum { value = true }; \ } AKANTU_SPECIFY_IS_SCALAR(Real); AKANTU_SPECIFY_IS_SCALAR(UInt); AKANTU_SPECIFY_IS_SCALAR(Int); AKANTU_SPECIFY_IS_SCALAR(bool); template < typename T1, typename T2 > struct is_same { enum { value = false }; // is_same represents a bool. }; template < typename T > struct is_same { enum { value = true }; }; /* -------------------------------------------------------------------------- */ #define AKANTU_SET_MACRO(name, variable, type) \ inline void set##name (type variable) { \ this->variable = variable; \ } #define AKANTU_GET_MACRO(name, variable, type) \ inline type get##name () const { \ return variable; \ } #define AKANTU_GET_MACRO_NOT_CONST(name, variable, type) \ inline type get##name () { \ return variable; \ } #define AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, \ support, con) \ inline con Array & \ get##name (const support & el_type, \ const GhostType & ghost_type = _not_ghost) con { \ return variable(el_type, ghost_type); \ } #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType,) #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, const) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType,) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, const) /* -------------------------------------------------------------------------- */ /// initialize the static part of akantu void initialize(int & argc, char ** & argv); /// initialize the static part of akantu and read the global input_file void initialize(const std::string & input_file, int & argc, char ** & argv); /* -------------------------------------------------------------------------- */ /// finilize correctly akantu and clean the memory void finalize (); /* -------------------------------------------------------------------------- */ /* * For intel compiler annoying remark */ #if defined(__INTEL_COMPILER) /// remark #981: operands are evaluated in unspecified order #pragma warning ( disable : 981 ) /// remark #383: value copied to temporary, reference to temporary used #pragma warning ( disable : 383 ) #endif //defined(__INTEL_COMPILER) /* -------------------------------------------------------------------------- */ /* string manipulation */ /* -------------------------------------------------------------------------- */ inline std::string to_lower(const std::string & str); /* -------------------------------------------------------------------------- */ inline std::string trim(const std::string & to_trim); /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /// give a string representation of the a human readable size in bit template std::string printMemorySize(UInt size); /* -------------------------------------------------------------------------- */ __END_AKANTU__ #include "aka_fwd.hh" __BEGIN_AKANTU__ /// get access to the internal argument parser cppargparse::ArgumentParser & getStaticArgumentParser(); /// get access to the internal input file parser const Parser & getStaticParser(); /// get access to the user part of the internal input file parser const ParserSection & getUserParser(); __END_AKANTU__ /* -------------------------------------------------------------------------- */ // BOOST PART: TOUCH ONLY IF YOU KNOW WHAT YOU ARE DOING #include #define AKANTU_BOOST_CASE_MACRO(r,macro,type) \ case type : { macro(type); break; } #define AKANTU_BOOST_ELEMENT_SWITCH(macro1, list1) \ do { \ switch(type) { \ BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_CASE_MACRO, macro1, list1) \ default: { \ AKANTU_DEBUG_ERROR("Type (" << UInt(type) << ") not handled by this function"); \ } \ } \ } while(0) #define AKANTU_BOOST_ALL_ELEMENT_SWITCH(macro) \ AKANTU_BOOST_ELEMENT_SWITCH(macro, \ AKANTU_ALL_ELEMENT_TYPE) #define AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(macro) \ AKANTU_BOOST_ELEMENT_SWITCH(macro, \ AKANTU_REGULAR_ELEMENT_TYPE) #define AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(macro) \ AKANTU_BOOST_ELEMENT_SWITCH(macro, \ AKANTU_COHESIVE_ELEMENT_TYPE) #define AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(macro) \ AKANTU_BOOST_ELEMENT_SWITCH(macro, \ AKANTU_STRUCTURAL_ELEMENT_TYPE) +#define AKANTU_BOOST_IGFEM_ELEMENT_SWITCH(macro) \ + AKANTU_BOOST_ELEMENT_SWITCH(macro, \ + AKANTU_IGFEM_ELEMENT_TYPE) + #define AKANTU_BOOST_LIST_MACRO(r,macro,type) \ macro(type) #define AKANTU_BOOST_ELEMENT_LIST(macro, list) \ BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_LIST_MACRO,macro,list) #define AKANTU_BOOST_ALL_ELEMENT_LIST(macro) \ AKANTU_BOOST_ELEMENT_LIST(macro, AKANTU_ALL_ELEMENT_TYPE) #define AKANTU_BOOST_REGULAR_ELEMENT_LIST(macro) \ AKANTU_BOOST_ELEMENT_LIST(macro, AKANTU_REGULAR_ELEMENT_TYPE) #define AKANTU_BOOST_STRUCTURAL_ELEMENT_LIST(macro) \ AKANTU_BOOST_ELEMENT_LIST(macro, AKANTU_STRUCTURAL_ELEMENT_TYPE) #define AKANTU_BOOST_COHESIVE_ELEMENT_LIST(macro) \ AKANTU_BOOST_ELEMENT_LIST(macro, AKANTU_COHESIVE_ELEMENT_TYPE) +#define AKANTU_BOOST_IGFEM_ELEMENT_LIST(macro) \ + AKANTU_BOOST_ELEMENT_LIST(macro, AKANTU_IGFEM_ELEMENT_TYPE) + #include "aka_common_inline_impl.cc" #endif /* __AKANTU_COMMON_HH__ */ diff --git a/src/common/aka_config.hh.in b/src/common/aka_config.hh.in index 843b4efc8..8d0613212 100644 --- a/src/common/aka_config.hh.in +++ b/src/common/aka_config.hh.in @@ -1,85 +1,87 @@ /** * @file aka_config.hh.in * * @author Nicolas Richart * * @date Fri Jan 13 12:34:54 2012 * * @brief Compilation time configuration of Akantu * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_AKA_CONFIG_HH__ #define __AKANTU_AKA_CONFIG_HH__ #define AKANTU_VERSION @AKANTU_VERSION@ #define AKANTU_VERSION_MAJOR @AKANTU_MAJOR_VERSION@ #define AKANTU_VERSION_MINOR @AKANTU_MINOR_VERSION@ #define AKANTU_VERSION_PATCH @AKANTU_PATCH_VERSION@ #cmakedefine AKANTU_USE_BLAS #cmakedefine AKANTU_USE_LAPACK #cmakedefine AKANTU_PARALLEL #cmakedefine AKANTU_USE_MPI #cmakedefine AKANTU_USE_SCOTCH #cmakedefine AKANTU_USE_PTSCOTCH #cmakedefine AKANTU_SCOTCH_NO_EXTERN #cmakedefine AKANTU_USE_MUMPS #cmakedefine AKANTU_USE_IOHELPER #cmakedefine AKANTU_USE_QVIEW #cmakedefine AKANTU_USE_BLACKDYNAMITE #cmakedefine AKANTU_USE_NLOPT #cmakedefine AKANTU_USE_CPPARRAY #cmakedefine AKANTU_USE_OBSOLETE_GETTIMEOFDAY #cmakedefine AKANTU_EXTRA_MATERIALS #cmakedefine AKANTU_DAMAGE_NON_LOCAL #cmakedefine AKANTU_STRUCTURAL_MECHANICS #cmakedefine AKANTU_HEAT_TRANSFER #cmakedefine AKANTU_COHESIVE_ELEMENT #cmakedefine AKANTU_PARALLEL_COHESIVE_ELEMENT +#cmakedefine AKANTU_IGFEM + // BOOST Section #cmakedefine AKANTU_BOOST_CHRONO #cmakedefine AKANTU_BOOST_SYSTEM // Experimental part #cmakedefine AKANTU_CORE_CXX11 // Debug tools #cmakedefine AKANTU_NDEBUG #cmakedefine AKANTU_DEBUG_TOOLS #cmakedefine READLINK_COMMAND @READLINK_COMMAND@ #cmakedefine ADDR2LINE_COMMAND @ADDR2LINE_COMMAND@ #define __aka_inline__ inline #endif /* __AKANTU_AKA_CONFIG_HH__ */ diff --git a/src/fem/element_class_igfem.hh b/src/fem/element_class_igfem.hh index 6a0750638..cdcdc080e 100644 --- a/src/fem/element_class_igfem.hh +++ b/src/fem/element_class_igfem.hh @@ -1,238 +1,270 @@ /** * @file element_class_igfem.hh * @author Aurelia Cuba Ramos * @author Nicolas Richart * @date Thu May 22 13:50:46 2014 * * @brief Specialization for interface-enriched finite elements * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ /* Interpolation */ /* -------------------------------------------------------------------------- */ template class InterpolationElement { public: typedef InterpolationPorperty interpolation_property; /// compute the shape values for a given set of points in natural coordinates static inline void computeShapes(const Matrix & natural_coord, Matrix & N) { - ElementClass::regular_el_type>::computeShapes(natural_coord, N); + AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute the shape values for a given point in natural coordinates template static inline void computeShapes(__attribute__((unused)) const vector_type & natural_coord, __attribute__((unused)) vector_type & N) { AKANTU_DEBUG_TO_IMPLEMENT(); } /** * compute @f$ B_{ij} = \frac{\partial N_j}{\partial S_i} @f$ the variation of * shape functions along with variation of natural coordinates on a given set * of points in natural coordinates */ static inline void computeDNDS(const Matrix & natural_coord, Tensor3 & dnds) { - ElementClass::regular_el_type>::computeShapes(natural_coord, N); + AKANTU_DEBUG_TO_IMPLEMENT(); } /** * compute @f$ B_{ij} = \frac{\partial N_j}{\partial S_i} @f$ the variation of shape functions along with * variation of natural coordinates on a given point in natural * coordinates */ template static inline void computeDNDS(__attribute__((unused)) const vector_type & natural_coord, __attribute__((unused)) matrix_type & dnds) { AKANTU_DEBUG_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(__attribute__((unused)) const Matrix & J, __attribute__((unused)) Real & jacobians) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// interpolate a field given (arbitrary) natural coordinates static inline void interpolateOnNaturalCoordinates(const Vector & natural_coords, const Matrix & nodal_values, Vector & interpolated) { - ElementClass::regular_el_type>::computeShapes(natural_coord, N); + AKANTU_DEBUG_TO_IMPLEMENT(); } /// 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) { - ElementClass::regular_el_type>::computeShapes(natural_coord, N); + AKANTU_DEBUG_TO_IMPLEMENT(); } public: static AKANTU_GET_MACRO_NOT_CONST(ShapeSize, InterpolationPorperty::nb_nodes_per_element, UInt); static AKANTU_GET_MACRO_NOT_CONST(ShapeDerivativesSize, (InterpolationPorperty::nb_nodes_per_element * InterpolationPorperty::natural_space_dimension), UInt); static AKANTU_GET_MACRO_NOT_CONST(NaturalSpaceDimension, InterpolationPorperty::natural_space_dimension, UInt); static AKANTU_GET_MACRO_NOT_CONST(NbNodesPerInterpolationElement, InterpolationPorperty::nb_nodes_per_element, UInt); }; /* -------------------------------------------------------------------------- */ /// define ElementClassProperty for parent elements #define AKANTU_DEFINE_IGFEM_PARENT_ELEMENT_CLASS_PROPERTY(elem_type, \ geom_type, \ interp_type, \ regular_el_type, \ - parent_el_type, \ elem_kind, \ sp, \ gauss_int_type, \ min_int_order) \ template<> \ struct ElementClassProperty { \ static const GeometricalType geometrical_type = geom_type; \ static const InterpolationType interpolation_type = interp_type; \ static const ElementType regular_element_type = regular_el_type; \ static const ElementKind element_kind = elem_kind; \ static const UInt spatial_dimension = sp; \ static const GaussIntergrationType gauss_integration_type = gauss_int_type; \ static const UInt minimal_integration_order = min_int_order; \ static const bool is_subelement = false; \ } /// define ElementClassProperty for subelements #define AKANTU_DEFINE_IGFEM_SUBELEMENT_CLASS_PROPERTY(elem_type, \ geom_type, \ interp_type, \ regular_el_type, \ parent_el_type, \ elem_kind, \ sp, \ gauss_int_type, \ min_int_order) \ template<> \ struct ElementClassProperty { \ static const GeometricalType geometrical_type = geom_type; \ static const InterpolationType interpolation_type = interp_type; \ static const ElementType regular_element_type = regular_el_type; \ static const ElementType parent_element_type = parent_el_type; \ static const ElementKind element_kind = elem_kind; \ static const UInt spatial_dimension = sp; \ static const GaussIntergrationType gauss_integration_type = gauss_int_type; \ static const UInt minimal_integration_order = min_int_order; \ static const bool is_subelement = true; \ } /* -------------------------------------------------------------------------- */ /* ElementClass */ /* -------------------------------------------------------------------------- */ template class ElementClass : public GeometricalElement::geometrical_type>, public InterpolationElement::interpolation_type> { protected: typedef GeometricalElement::geometrical_type> geometrical_element; typedef InterpolationElement::interpolation_type> interpolation_element; typedef ElementClassProperty element_property; typedef typename interpolation_element::interpolation_property interpolation_property; public: + /// compute the shape values for a given set of points in natural coordinates + static inline void computeShapes(const Matrix & natural_coord, + Matrix & N) { + ElementClass::regular_el_type>::computeShapes(natural_coord, N); + } + + static inline void computeDNDS(const Matrix & natural_coord, + Tensor3 & dnds) { + ElementClass::regular_el_type>::computeDNDS(natural_coord, dnds); + } + + /// interpolate a field given (arbitrary) natural coordinates + static inline void interpolateOnNaturalCoordinates(const Vector & natural_coords, + const Matrix & nodal_values, + Vector & interpolated) { + ElementClass::regular_el_type>::interpolateonnaturalcoordinates(natural_coords, nodal_values, 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) { + ElementClass::regular_el_type>::gradientOnNaturalCoordinates(natural_coords, f, gradient); + } + + /** * 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) { - ElementClass::regular_el_type>::computeShapes(natural_coord, N); + ElementClass::regular_el_type>::computeJMat(dnds, node_coords, 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) { - ElementClass::regular_el_type>::computeShapes(natural_coord, N); + ElementClass::regular_el_type>::JMat(dnds, node_coords, J); } /// compute the jacobians of a serie of natural coordinates static inline void computeJacobian(const Matrix & natural_coords, const Matrix & node_coords, Vector & jacobians) { - ElementClass::regular_el_type>::computeShapes(natural_coord, N); + ElementClass::regular_el_type>::computeJacobian(natural_coords, node_coords, jacobians); } /// compute jacobian (or integration variable change factor) for a set of points static inline void computeJacobian(const Tensor3 & J, Vector & jacobians) { - ElementClass::regular_el_type>::computeShapes(natural_coord, N); + ElementClass::regular_el_type>::computeJacobian(J, jacobians); } /// compute jacobian (or integration variable change factor) for a given point static inline void computeJacobian(const Matrix & J, Real & jacobians) { - ElementClass::regular_el_type>::computeShapes(natural_coord, N); + ElementClass::regular_el_type>::computeJacobian(J, 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) { - ElementClass::regular_el_type>::computeShapes(natural_coord, N); + ElementClass::regular_el_type>::computeShapeDerivatives(J,dnds, 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); + Matrix & normals){ + ElementClass::regular_el_type>::computeNormalsOnNaturalCoordinates(coord, f, normals); + } /// get natural coordinates from real coordinates static inline void inverseMap(const Vector & real_coords, const Matrix & node_coords, Vector & natural_coords, Real tolerance = 1e-8) { - ElementClass::regular_el_type>::computeShapes(natural_coord, N); + ElementClass::regular_el_type>::inverseMap(real_coords, node_coords, natural_coords, tolerance); } public: - static AKANTU_GET_MACRO_NOT_CONST(Kind, element_kind, ElementKind); + static AKANTU_GET_MACRO_NOT_CONST(Kind, _ek_igfem, ElementKind); static AKANTU_GET_MACRO_NOT_CONST(SpatialDimension, ElementClassProperty::spatial_dimension, UInt); static AKANTU_GET_MACRO_NOT_CONST(P1ElementType, p1_type, const ElementType &); static const ElementType & getFacetType(UInt t = 0) { return facet_type[t]; } private: /// Type of the facet elements static ElementType facet_type[]; /// type of element P1 associated static ElementType p1_type; }; /* -------------------------------------------------------------------------- */ diff --git a/src/fem/igfem_element.hh b/src/fem/igfem_element.hh new file mode 100644 index 000000000..dba85e43b --- /dev/null +++ b/src/fem/igfem_element.hh @@ -0,0 +1,44 @@ +/** + * @file igfem_element.hh + * @author Aurelia Cuba Ramos + * @date Fri May 30 17:56:54 2014 + * + * @brief IGFEM element class + * + * @section LICENSE + * + * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) + * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) + * + * Akantu is free software: you can redistribute it and/or modify it under the + * terms of the GNU Lesser General Public License as published by the Free + * Software Foundation, either version 3 of the License, or (at your option) any + * later version. + * + * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR + * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more + * details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with Akantu. If not, see . + * + */ + +/* -------------------------------------------------------------------------- */ + +/* -------------------------------------------------------------------------- */ +#include "element_class_igfem.hh" + +/* -------------------------------------------------------------------------- */ + +#ifndef __AKANTU_IGFEM_ELEMENT_HH__ +#define __AKANTU_IGFEM_ELEMENT_HH__ + +__BEGIN_AKANTU__ + +AKANTU_DEFINE_IGFEM_PARENT_ELEMENT_CLASS_PROPERTY(_igfem_triangle_3, _gt_triangle_3, _itp_lagrange_triangle_3, _triangle_3, _ek_igfem, 2, + _git_triangle, 1); +__END_AKANTU__ + +#endif /* __AKANTU_IGFEM_ELEMENT_HH__ */ diff --git a/src/fem/shape_igfem.hh b/src/fem/shape_igfem.hh index 59b5c74aa..a2c0b3b85 100644 --- a/src/fem/shape_igfem.hh +++ b/src/fem/shape_igfem.hh @@ -1,201 +1,202 @@ /** * @file shape_igfem.hh * @author Aurelia Cuba Ramos * @date Sun May 25 10:34:41 2014 * * @brief shape functions for interface-enriched generalized FEM * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ #ifndef __AKANTU_SHAPE_IGFEM_HH__ #define __AKANTU_SHAPE_IGFEM_HH__ #include "shape_functions.hh" +#include "igfem_element.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template<> class ShapeLagrange<_ek_igfem> : public ShapeFunctions{ /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: ShapeLagrange(const Mesh & mesh, const ID & id = "shape_igfem", const MemoryID & memory_id = 0); virtual ~ShapeLagrange(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: inline void initShapeFunctions(const Array & nodes, const Matrix & control_points, const ElementType & type, const GhostType & ghost_type); /// set the control points for a given element template void setControlPointsByType(const Matrix & control_points, const GhostType & ghost_type); /// pre compute all shapes on the element control points from natural coordinates template void precomputeShapesOnControlPoints(const Array & nodes, GhostType ghost_type); /// pre compute all shape derivatives on the element control points from natural coordinates template void precomputeShapeDerivativesOnControlPoints(const Array & nodes, GhostType ghost_type); /// interpolate nodal values on the control points template void interpolateOnControlPoints(const Array &u, Array &uq, UInt nb_degree_of_freedom, GhostType ghost_type = _not_ghost, const Array & filter_elements = empty_filter) const; template void interpolateOnControlPoints(const Array &u, Array &uq, UInt nb_degree_of_freedom, GhostType ghost_type = _not_ghost, const Array & filter_elements = empty_filter) const; /// compute the gradient of u on the control points template void gradientOnControlPoints(const Array &u, Array &nablauq, UInt nb_degree_of_freedom, GhostType ghost_type = _not_ghost, const Array & filter_elements = empty_filter) const; template void gradientOnControlPoints(const Array &u, Array &nablauq, UInt nb_degree_of_freedom, GhostType ghost_type = _not_ghost, const Array & filter_elements = empty_filter) const; /// multiply a field by shape functions @f$ fts_{ij} = f_i * \varphi_j @f$ template void fieldTimesShapes(const Array & field, Array & field_times_shapes, GhostType ghost_type) const; template void fieldTimesShapes(const Array & field, Array & field_times_shapes, GhostType ghost_type) const; /// find natural coords from real coords provided an element template void inverseMap(const Vector & real_coords, UInt element, Vector & natural_coords, const GhostType & ghost_type = _not_ghost) const; template void inverseMap(const Vector & real_coords, UInt element, Vector & natural_coords, const GhostType & ghost_type = _not_ghost) const; /// return true if the coordinates provided are inside the element, false otherwise template bool contains(const Vector & real_coords, UInt elem, const GhostType & ghost_type) const; template bool contains(const Vector & real_coords, UInt elem, const GhostType & ghost_type) const; /// compute the shape on a provided point template void computeShapes(const Vector & real_coords, UInt elem, Vector & shapes, const GhostType & ghost_type) const; template void computeShapes(const Vector & real_coords, UInt elem, Vector & shapes, const GhostType & ghost_type) const; /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; template void printself(std::ostream & stream, int indent = 0) const; protected: /// compute the shape derivatives on control points for a given element template inline void computeShapeDerivativesOnCPointsByElement(const Matrix & node_coords, const Matrix & natural_coords, Tensor3 & shapesd); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get a the shapes vector inline const Array & getShapes(const ElementType & el_type, const GhostType & ghost_type = _not_ghost) const; /// get a the shapes derivatives vector inline const Array & getShapesDerivatives(const ElementType & el_type, const GhostType & ghost_type = _not_ghost) const; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// shape functions for all elements ByElementTypeArray shapes; /// shape functions derivatives for all elements ByElementTypeArray shapes_derivatives; /// for parent interpolation types store boolean to decide wheter element is enriched or not ByElementTypeArray is_split; /// for subelement store the control points of the corresponding parent elements ByElementTypeArray parent_control_points; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "shape_igfem_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_SHAPE_IGFEM_HH__ */ diff --git a/src/fem/shape_igfem_inline_impl.cc b/src/fem/shape_igfem_inline_impl.cc index 6ec2f7a5b..886103968 100644 --- a/src/fem/shape_igfem_inline_impl.cc +++ b/src/fem/shape_igfem_inline_impl.cc @@ -1,566 +1,551 @@ /** * @file shape_igfem_inline_impl.cc * @author Aurelia Cuba Ramos * @date Sun May 25 21:53:21 2014 * * @brief ShapeIGFEM inline implementation * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ -__END_AKANTU__ - -#include "fem.hh" - -__BEGIN_AKANTU__ - /* -------------------------------------------------------------------------- */ ShapeLagrange<_ek_igfem>::ShapeLagrange(const Mesh & mesh, const ID & id, const MemoryID & memory_id) : ShapeFunctions(mesh, id, memory_id), shapes("shapes_generic", id), shapes_derivatives("shapes_derivatives_generic", id) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -inline const Array & ShapeIGFEM::getShapes(const ElementType & el_type, +inline const Array & ShapeLagrange<_ek_igfem>::getShapes(const ElementType & el_type, const GhostType & ghost_type) const { return shapes(FEM::getInterpolationType(el_type), ghost_type); } /* -------------------------------------------------------------------------- */ inline const Array & ShapeLagrange<_ek_igfem>::getShapesDerivatives(const ElementType & el_type, const GhostType & ghost_type) const { return shapes_derivatives(FEM::getInterpolationType(el_type), ghost_type); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_igfem>::setControlPointsByType(const Matrix & points, const GhostType & ghost_type) { if (!is_sub) - ShapeFunctions<_ek_igfem>::control_points(type, ghost_type).shallowCopy(points); + ShapeLagrange<_ek_igfem>::control_points(type, ghost_type).shallowCopy(points); else { AKANTU_DEBUG_TO_IMPLEMENT(); } } /* -------------------------------------------------------------------------- */ #define INIT_SHAPE_FUNCTIONS(type) \ setControlPointsByType::is_subelement>(control_points, ghost_type); \ precomputeShapesOnControlPoints::is_subelement>(nodes, ghost_type); \ if (ElementClass::getNaturalSpaceDimension() == \ mesh.getSpatialDimension() || kind != _ek_regular) \ precomputeShapeDerivativesOnControlPoints::is_subelement>(nodes, ghost_type); -template -inline void -ShapeLagrange::initShapeFunctions(const Array & nodes, - const Matrix & control_points, - const ElementType & type, - const GhostType & ghost_type) { - AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(INIT_SHAPE_FUNCTIONS); -} - /* -------------------------------------------------------------------------- */ template inline void ShapeLagrange<_ek_igfem>:: computeShapeDerivativesOnCPointsByElement(const Matrix & node_coords, const Matrix & natural_coords, Tensor3 & shapesd) { if(!is_sub) { // // compute dnds Tensor3 dnds(node_coords.rows(), node_coords.cols(), natural_coords.cols()); ElementClass::computeDNDS(natural_coords, dnds); // compute dxds Tensor3 J(node_coords.rows(), natural_coords.rows(), natural_coords.cols()); ElementClass::computeJMat(dnds, node_coords, J); // compute shape derivatives ElementClass::computeShapeDerivatives(J, dnds, shapesd); } else { - AKANTU_DEBUG_TO_IMPLEMENT; + AKANTU_DEBUG_TO_IMPLEMENT(); } } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_igfem>::inverseMap(const Vector & real_coords, UInt elem, Vector & natural_coords, const GhostType & ghost_type) const{ - inverseMap::is_subelement>(real_coords, elem, natural_coods, ghost_type); + inverseMap::is_subelement>(real_coords, elem, natural_coords, ghost_type); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_igfem>::inverseMap(const Vector & real_coords, UInt elem, Vector & natural_coords, const GhostType & ghost_type) const{ if (!is_sub) { UInt spatial_dimension = mesh.getSpatialDimension(); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); UInt * elem_val = mesh.getConnectivity(type, ghost_type).storage(); Matrix nodes_coord(spatial_dimension, nb_nodes_per_element); mesh.extractNodalValuesFromElement(mesh.getNodes(), nodes_coord.storage(), elem_val + elem*nb_nodes_per_element, nb_nodes_per_element, spatial_dimension); ElementClass::inverseMap(real_coords, nodes_coord, natural_coords); } else { - AKANTU_DEBUG_TO_IMPLEMENT; + AKANTU_DEBUG_TO_IMPLEMENT(); } } /* -------------------------------------------------------------------------- */ template bool ShapeLagrange<_ek_igfem>::contains(const Vector & real_coords, UInt elem, const GhostType & ghost_type) const{ - return contains::is_subelement>(real_coods, elem, ghost_type); + return contains::is_subelement>(real_coords, elem, ghost_type); } /* -------------------------------------------------------------------------- */ template bool ShapeLagrange<_ek_igfem>::contains(const Vector & real_coords, UInt elem, const GhostType & ghost_type) const{ if (!is_sub) { UInt spatial_dimension = mesh.getSpatialDimension(); Vector natural_coords(spatial_dimension); inverseMap(real_coords, elem, natural_coords, ghost_type); return ElementClass::contains(natural_coords); } else { - AKANTU_DEBUG_TO_IMPLEMENT; + AKANTU_DEBUG_TO_IMPLEMENT(); } } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_igfem>::computeShapes(const Vector & real_coords, UInt elem, Vector & shapes, const GhostType & ghost_type) const{ computeShapes::is_subelement>(real_coords, elem, shapes, ghost_type); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_igfem>::computeShapes(const Vector & real_coords, UInt elem, Vector & shapes, const GhostType & ghost_type) const{ if (!is_sub) { UInt spatial_dimension = mesh.getSpatialDimension(); Vector natural_coords(spatial_dimension); inverseMap(real_coords, elem, natural_coords, ghost_type); ElementClass::computeShapes(natural_coords, shapes); } else { - AKANTU_DEBUG_TO_IMPLEMENT; + AKANTU_DEBUG_TO_IMPLEMENT(); } } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_igfem>::precomputeShapesOnControlPoints(__attribute__((unused)) const Array & nodes, GhostType ghost_type) { AKANTU_DEBUG_IN(); if (!is_sub) { InterpolationType itp_type = ElementClassProperty::interpolation_type; Matrix & natural_coords = control_points(type, ghost_type); UInt nb_points = natural_coords.cols(); UInt size_of_shapes = ElementClass::getShapeSize(); UInt nb_element = mesh.getConnectivity(type, ghost_type).getSize();; Array & shapes_tmp = shapes.alloc(nb_element*nb_points, size_of_shapes, itp_type, ghost_type); Array::matrix_iterator shapes_it = shapes_tmp.begin_reinterpret(ElementClass::getNbNodesPerInterpolationElement(), nb_points, nb_element); for (UInt elem = 0; elem < nb_element; ++elem, ++shapes_it) { Matrix & N = *shapes_it; ElementClass::computeShapes(natural_coords, N); } } else { - AKANTU_DEBUG_TO_IMPLEMENT; + AKANTU_DEBUG_TO_IMPLEMENT(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_igfem>::precomputeShapeDerivativesOnControlPoints(const Array & nodes, GhostType ghost_type) { AKANTU_DEBUG_IN(); if (!is_sub) { InterpolationType itp_type = ElementClassProperty::interpolation_type; // Real * coord = mesh.getNodes().storage(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); UInt size_of_shapesd = ElementClass::getShapeDerivativesSize(); Matrix & natural_coords = control_points(type, ghost_type); UInt nb_points = natural_coords.cols(); UInt nb_element = mesh.getConnectivity(type, ghost_type).getSize(); Array & shapes_derivatives_tmp = shapes_derivatives.alloc(nb_element*nb_points, size_of_shapesd, itp_type, ghost_type); Array x_el(0, spatial_dimension * nb_nodes_per_element); FEM::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type); Real * shapesd_val = shapes_derivatives_tmp.storage(); Array::matrix_iterator x_it = x_el.begin(spatial_dimension, nb_nodes_per_element); for (UInt elem = 0; elem < nb_element; ++elem, ++x_it) { Matrix & X = *x_it; Tensor3 B(shapesd_val, spatial_dimension, nb_nodes_per_element, nb_points); computeShapeDerivativesOnCPointsByElement::is_subelement>(X, natural_coords, B); shapesd_val += size_of_shapesd*nb_points; } } else { - AKANTU_DEBUG_TO_IMPLEMENT; + AKANTU_DEBUG_TO_IMPLEMENT(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_igfem>::interpolateOnControlPoints(const Array &in_u, Array &out_uq, UInt nb_degree_of_freedom, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); interpolateOnControlPoints::is_subelement>(in_u, out_uq, nb_degree_of_freedom, ghost_type, filter_elements); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_igfem>::interpolateOnControlPoints(const Array &in_u, Array &out_uq, UInt nb_degree_of_freedom, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); if(!is_sub) { InterpolationType itp_type = ElementClassProperty::interpolation_type; AKANTU_DEBUG_ASSERT(shapes.exists(itp_type, ghost_type), "No shapes for the type " << shapes.printType(itp_type, ghost_type)); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); Array u_el(0, nb_degree_of_freedom * nb_nodes_per_element); FEM::extractNodalToElementField(mesh, in_u, u_el, type, ghost_type, filter_elements); this->interpolateElementalFieldOnControlPoints(u_el, out_uq, ghost_type, shapes(itp_type, ghost_type), filter_elements); } else { - AKANTU_DEBUG_TO_IMPLEMENT; + AKANTU_DEBUG_TO_IMPLEMENT(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_igfem>::gradientOnControlPoints(const Array &in_u, Array &out_nablauq, UInt nb_degree_of_freedom, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); gradientOnControlPoints::is_subelement>(in_u, out_nablauq, nb_degree_of_freedom, ghost_type, filter_elements); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_igfem>::gradientOnControlPoints(const Array &in_u, Array &out_nablauq, UInt nb_degree_of_freedom, GhostType ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); if (!is_sub) { InterpolationType itp_type = ElementClassProperty::interpolation_type; AKANTU_DEBUG_ASSERT(shapes_derivatives.exists(itp_type, ghost_type), "No shapes derivatives for the type " << shapes_derivatives.printType(itp_type, ghost_type)); UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); Array u_el(0, nb_degree_of_freedom * nb_nodes_per_element); FEM::extractNodalToElementField(mesh, in_u, u_el, type, ghost_type, filter_elements); this->gradientElementalFieldOnControlPoints(u_el, out_nablauq, ghost_type, shapes_derivatives(itp_type, ghost_type), filter_elements); } else { - AKANTU_DEBUG_ASSERT; + AKANTU_DEBUG_TO_IMPLEMENT(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_igfem>::fieldTimesShapes(const Array & field, Array & field_times_shapes, GhostType ghost_type) const { fieldTimesShapes::is_subelement>(field, field_times_shapes, ghost_type); } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_igfem>::fieldTimesShapes(const Array & field, Array & field_times_shapes, GhostType ghost_type) const { if (!is_sub) { field_times_shapes.resize(field.getSize()); UInt size_of_shapes = ElementClass::getShapeSize(); InterpolationType itp_type = ElementClassProperty::interpolation_type; UInt nb_degree_of_freedom = field.getNbComponent(); const Array & shape = shapes(itp_type, ghost_type); Array::const_matrix_iterator field_it = field.begin(nb_degree_of_freedom, 1); Array::const_matrix_iterator shapes_it = shape.begin(1, size_of_shapes); Array::matrix_iterator it = field_times_shapes.begin(nb_degree_of_freedom, size_of_shapes); Array::matrix_iterator end = field_times_shapes.end (nb_degree_of_freedom, size_of_shapes); for (; it != end; ++it, ++field_it, ++shapes_it) { it->mul(*field_it, *shapes_it); } } else { - AKANTU_DEBUG_TO_IMPLEMENT; + AKANTU_DEBUG_TO_IMPLEMENT(); } } /* -------------------------------------------------------------------------- */ void ShapeLagrange<_ek_igfem>::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "Shapes Lagrange [" << std::endl; ShapeFunctions::printself(stream, indent + 1); shapes.printself(stream, indent + 1); shapes_derivatives.printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ template void ShapeLagrange<_ek_igfem>::printself(std::ostream & stream, int indent) const { if (!is_sub) { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "Shapes Lagrange [" << std::endl; ShapeFunctions::printself(stream, indent + 1); shapes.printself(stream, indent + 1); shapes_derivatives.printself(stream, indent + 1); stream << space << "]" << std::endl; } else { - AKANTU_DEBUG_TO_IMPLEMENT; + AKANTU_DEBUG_TO_IMPLEMENT(); } } /* -------------------------------------------------------------------------- */ // __END_AKANTU__ // #include "fem.hh" // __BEGIN_AKANTU__ // /* -------------------------------------------------------------------------- */ // template // template // void ShapeLagrange::precomputeShapesOnControlPoints(__attribute__((unused)) const Array & nodes, // GhostType ghost_type) { // AKANTU_DEBUG_IN(); // InterpolationType itp_type = ElementClassProperty::interpolation_type; // /// get natural coordinates of subelement // Matrix & natural_coords = control_points(type, ghost_type); // UInt nb_points = natural_coords.cols(); // UInt size_of_shapes = ElementClass::getShapeSize(); // UInt nb_element = mesh.getConnectivity(type, ghost_type).getSize();; // Array & shapes_tmp = shapes.alloc(nb_element*nb_points, // size_of_shapes, // itp_type, // ghost_type); // Array::matrix_iterator shapes_it = // shapes_tmp.begin_reinterpret(ElementClass::getNbNodesPerInterpolationElement(), nb_points, // nb_element); // Array & parent_control_points_tmp = parent_control_points.alloc(nb_element, // nb_points, // itp_type, // ghost_type); // Array::matrix_iterator parent_control_points_it = // parent_control_points_tmp.begin(nb_points,spatial_dimension); // Matrix physical_control_points(nb_points,nb_points); // physical_control_points.clear(); // Matrix N_tmp = (); // for (UInt elem = 0; elem < nb_element; ++elem, ++shapes_it; ++parent_control_points_it) { // Matrix & N = *shapes_it; // ElementClass::computeShapes(natural_coords, // N_tmp); // if (!ElementClassProperty::has_parent) // N = N_tmp; // else { // /// compute control points in physical domain // physical_control_points.mul(natural_coords, control_points); // Array::vector_iterator physical_control_points_it = physical_control_points.begin(spatial_dimension); // /// compute control points in reference domain of the parent // for (UInt i = 0; i < spatial_dimension; ++i) { // inverseMap(physical_control_points, elem, *physical_control_points_it, ghost_type); // } // /// compute parent shape functions // ElementClass::parent_el_type>::computeShapes(natural_coords_parents) // } // } // AKANTU_DEBUG_OUT(); // } // /* -------------------------------------------------------------------------- */ // template // template // void ShapeLagrange::inverseMap(const Vector & real_coords, // UInt elem, // Vector & natural_coords, // const GhostType & ghost_type) const{ // UInt spatial_dimension = mesh.getSpatialDimension(); // UInt nb_nodes_per_element = ElementClass::getNbNodesPerInterpolationElement(); // UInt * elem_val = mesh.getConnectivity(type, ghost_type).storage(); // Matrix nodes_coord(spatial_dimension, nb_nodes_per_element); // mesh.extractNodalValuesFromElement(mesh.getNodes(), // nodes_coord.storage(), // elem_val + elem*nb_nodes_per_element, // nb_nodes_per_element, // spatial_dimension); // ElementClass::inverseMap(real_coords, // nodes_coord, // natural_coords); // } // /* -------------------------------------------------------------------------- */ diff --git a/test/test_fem/test_igfem_integrate.cc b/test/test_fem/test_igfem_integrate.cc new file mode 100644 index 000000000..5b1c75326 --- /dev/null +++ b/test/test_fem/test_igfem_integrate.cc @@ -0,0 +1,75 @@ +#include "aka_common.hh" +#include "shape_igfem.hh" +#include "integrator_gauss.hh" +/* -------------------------------------------------------------------------- */ +#include +#include +#include +/* -------------------------------------------------------------------------- */ + +using namespace akantu; + +int main(int argc, char *argv[]) { + akantu::initialize(argc, argv); + debug::setDebugLevel(dblTest); + const ElementType type = _igfem_triangle_3; + UInt dim = ElementClass::getSpatialDimension(); + + Real eps = 3e-13; + std::cout << "Epsilon : " << eps << std::endl; + + MeshIOMSH mesh_io; + Mesh my_mesh(dim); + + std::stringstream meshfilename; meshfilename << "_triangle_3.msh"; + mesh_io.read(meshfilename.str(), my_mesh); + + FEM *fem = new FEMTemplate >(my_mesh, dim, "my_fem"); + + std::stringstream outfilename; outfilename << "out_" << type << ".txt"; + std::ofstream my_file(outfilename.str().c_str()); + + fem->initShapeFunctions(); + + std::cout << *fem << std::endl; + + UInt nb_element = my_mesh.getNbElement(type); + UInt nb_quadrature_points = fem->getNbQuadraturePoints(type) * nb_element; + + Array const_val(fem->getMesh().getNbNodes(), 2, "const_val"); + Array val_on_quad(nb_quadrature_points, 2 , "val_on_quad"); + + for (UInt i = 0; i < const_val.getSize(); ++i) { + const_val.storage()[i * 2 + 0] = 1.; + const_val.storage()[i * 2 + 1] = 2.; + } + + //interpolate function on quadrature points + fem->interpolateOnQuadraturePoints(const_val, val_on_quad, 2, type); + + //integrate function on elements + akantu::Array int_val_on_elem(nb_element, 2, "int_val_on_elem"); + fem->integrate(val_on_quad, int_val_on_elem, 2, type); + + // get global integration value + Real value[2] = {0,0}; + my_file << val_on_quad << std::endl << int_val_on_elem << std::endl; + for (UInt i = 0; i < fem->getMesh().getNbElement(type); ++i) { + value[0] += int_val_on_elem.storage()[2*i]; + value[1] += int_val_on_elem.storage()[2*i+1]; + } + + my_file << "integral on the mesh of 1 is " << value[0] << " and of 2 is " << value[1] << std::endl; + + + delete fem; + finalize(); + + if(!(std::abs(value[0] - 1.) < eps && std::abs(value[1] - 2.) < eps)) { + std::cout << "|1 - " << value[0] << "| = " << std::abs(value[0] - 1.) << std::endl + << "|2 - " << value[1] << "| = " << std::abs(value[1] - 2.) << std::endl; + return EXIT_FAILURE; + } + + return EXIT_SUCCESS; +}