diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh index dde1e5e35..a66f1a43f 100644 --- a/src/common/aka_common.hh +++ b/src/common/aka_common.hh @@ -1,494 +1,495 @@ /** * @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 UInt Surface; typedef std::pair SurfacePair; typedef std::list< SurfacePair > SurfacePairList; /* -------------------------------------------------------------------------- */ /// @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) \ (_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) #else # define AKANTU_COHESIVE_ELEMENT_TYPE #endif #define AKANTU_ALL_ELEMENT_TYPE \ AKANTU_REGULAR_ELEMENT_TYPE \ AKANTU_COHESIVE_ELEMENT_TYPE \ AKANTU_STRUCTURAL_ELEMENT_TYPE #define AKANTU_NOT_STRUCTURAL_ELEMENT_TYPE \ AKANTU_REGULAR_ELEMENT_TYPE \ AKANTU_COHESIVE_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 _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 #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 #if defined(AKANTU_COHESIVE_ELEMENT) _gt_cohesive_2d_4, ///< 4 nodes 2D cohesive _gt_cohesive_2d_6, ///< 6 nodes 2D cohesive #endif _gt_not_defined }; /// @enum InterpolationType type of elements enum InterpolationType { _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_bernoulli_beam, ///< Bernoulli beam _itp_not_defined }; /// @enum InterpolationKind the familly of interpolation types enum InterpolationKind { _itk_lagrangian, _itk_structural }; //! standard output stream operator for ElementType inline std::ostream & operator <<(std::ostream & stream, ElementType type); enum ElementKind { _ek_regular, _ek_cohesive, _ek_structural, _ek_not_defined }; /// enum MeshIOType type of mesh reader/writer enum MeshIOType { _miot_auto, _miot_gmsh, _miot_diana }; /// enum AnalysisMethod type of solving method used to solve the equation of motion enum AnalysisMethod { _static, _implicit_dynamic, _explicit_dynamic }; /// 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, //< synchronization of facet data to setup facet synch + _gst_smmc_normals, //< synchronization of facet normals to setup facet synch //--- 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 }; /// 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 /* -------------------------------------------------------------------------- */ 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) /* -------------------------------------------------------------------------- */ void initialize(int & argc, char ** & argv); /* -------------------------------------------------------------------------- */ 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); __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 (" << 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_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) #include "aka_common_inline_impl.cc" #include "aka_fwd.hh" #endif /* __AKANTU_COMMON_HH__ */ diff --git a/src/common/aka_common_inline_impl.cc b/src/common/aka_common_inline_impl.cc index fd4a7cb9f..53d007356 100644 --- a/src/common/aka_common_inline_impl.cc +++ b/src/common/aka_common_inline_impl.cc @@ -1,136 +1,137 @@ /** * @file aka_common_inline_impl.cc * * @author Guillaume Anciaux * @author Nicolas Richart * * @date Thu Dec 01 12:54:29 2011 * * @brief inline implementations of common akantu type descriptions * * @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 * */ #include __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ //! standard output stream operator for ElementType inline std::ostream & operator <<(std::ostream & stream, ElementType type) { #define STRINGIFY(type) \ stream << BOOST_PP_STRINGIZE(type); AKANTU_BOOST_ALL_ELEMENT_SWITCH(STRINGIFY); #undef STRINGIFY return stream; } /// standard output stream operator for InterpolationType inline std::ostream & operator <<(std::ostream & stream, InterpolationType type) { switch(type) { case _itp_lagrange_segment_2 : stream << "_itp_lagrange_segment_2" ; break; case _itp_lagrange_segment_3 : stream << "_itp_lagrange_segment_3" ; break; case _itp_lagrange_triangle_3 : stream << "_itp_lagrange_triangle_3" ; break; case _itp_lagrange_triangle_6 : stream << "_itp_lagrange_triangle_6" ; break; case _itp_lagrange_quadrangle_4 : stream << "_itp_lagrange_quadrangle_4" ; break; case _itp_serendip_quadrangle_8 : stream << "_itp_serendip_quadrangle_8" ; break; case _itp_lagrange_tetrahedron_4 : stream << "_itp_lagrange_tetrahedron_4" ; break; case _itp_lagrange_tetrahedron_10 : stream << "_itp_lagrange_tetrahedron_10"; break; case _itp_lagrange_hexahedron_8 : stream << "_itp_lagrange_hexahedron_8" ; break; case _itp_bernoulli_beam : stream << "_itp_bernoulli_beam" ; break; case _itp_not_defined : stream << "_itp_not_defined" ; break; } return stream; } /// standard output stream operator for GhostType inline std::ostream & operator <<(std::ostream & stream, GhostType type) { switch(type) { case _not_ghost : stream << "not_ghost"; break; case _ghost : stream << "ghost" ; break; case _casper : stream << "Casper the friendly ghost"; break; } return stream; } /// standard output stream operator for SynchronizationTag inline std::ostream & operator <<(std::ostream & stream, SynchronizationTag type) { switch(type) { case _gst_smm_mass : stream << "_gst_smm_mass" ; break; case _gst_smm_for_strain : stream << "_gst_smm_for_strain" ; break; case _gst_smm_boundary : stream << "_gst_smm_boundary" ; break; case _gst_smm_uv : stream << "_gst_smm_uv" ; break; case _gst_smm_res : stream << "_gst_smm_res" ; break; case _gst_smm_init_mat : stream << "_gst_smm_init_mat" ; break; case _gst_smm_stress : stream << "_gst_smm_stress" ; break; case _gst_smmc_facets : stream << "_gst_smmc_facets" ; break; + case _gst_smmc_normals : stream << "_gst_smmc_normals" ; break; case _gst_htm_capacity : stream << "_gst_htm_capacity" ; break; case _gst_htm_temperature : stream << "_gst_htm_temperature" ; break; case _gst_htm_gradient_temperature : stream << "_gst_htm_gradient_temperature"; break; case _gst_htm_phi : stream << "_gst_htm_phi" ; break; case _gst_htm_gradient_phi : stream << "_gst_htm_gradient_phi" ; break; case _gst_mnl_for_average : stream << "_gst_mnl_for_average" ; break; case _gst_mnl_weight : stream << "_gst_mnl_weight" ; break; case _gst_test : stream << "_gst_test" ; break; case _gst_material_id : stream << "_gst_material_id" ; break; } return stream; } /* -------------------------------------------------------------------------- */ inline std::string to_lower(const std::string & str) { std::string lstr = str; std::transform(lstr.begin(), lstr.end(), lstr.begin(), (int(*)(int))std::tolower); return lstr; } /* -------------------------------------------------------------------------- */ inline std::string trim(const std::string & to_trim) { std::string trimed = to_trim; //left trim trimed.erase(trimed.begin(), std::find_if(trimed.begin(), trimed.end(), std::not1(std::ptr_fun(std::isspace)))); // right trim trimed.erase(std::find_if(trimed.rbegin(), trimed.rend(), std::not1(std::ptr_fun(std::isspace))).base(), trimed.end()); return trimed; } __END_AKANTU__ diff --git a/src/common/aka_vector.hh b/src/common/aka_vector.hh index 55f2942a0..dc053d93e 100644 --- a/src/common/aka_vector.hh +++ b/src/common/aka_vector.hh @@ -1,369 +1,372 @@ /** * @file aka_vector.hh * * @author Nicolas Richart * * @date Fri Jun 18 11:48:28 2010 * * @brief class of vectors * * @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_VECTOR_HH__ #define __AKANTU_VECTOR_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" /* -------------------------------------------------------------------------- */ #include //#include #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /// class that afford to store vectors in static memory class ArrayBase { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: ArrayBase(const ID & id = ""); virtual ~ArrayBase(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// get the amount of space allocated in bytes inline UInt getMemorySize() const; /// set the size to zero without freeing the allocated space inline void empty(); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(AllocatedSize, allocated_size, UInt); AKANTU_GET_MACRO(Size, size, UInt); AKANTU_GET_MACRO(NbComponent, nb_component, UInt); AKANTU_GET_MACRO(ID, id, const ID &); AKANTU_GET_MACRO(Tag, tag, const std::string &); AKANTU_SET_MACRO(Tag, tag, const std::string &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id of the vector ID id; /// the size allocated UInt allocated_size; /// the size used UInt size; /// number of components UInt nb_component; /// size of the stored type UInt size_of_type; /// User defined tag std::string tag; }; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template class Array : public ArrayBase { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef T value_type; typedef value_type & reference; typedef value_type * pointer_type; typedef const value_type & const_reference; /// Allocation of a new vector inline Array(UInt size = 0, UInt nb_component = 1, const ID & id = ""); /// Allocation of a new vector with a default value Array(UInt size, UInt nb_component, const value_type def_values[], const ID & id = ""); /// Allocation of a new vector with a default value Array(UInt size, UInt nb_component, const_reference value, const ID & id = ""); /// Copy constructor (deep copy if deep=true) Array(const Array& vect, bool deep = true, const ID & id = ""); /// Copy constructor (deep copy) Array(const std::vector & vect); virtual inline ~Array(); /* ------------------------------------------------------------------------ */ /* Iterator */ /* ------------------------------------------------------------------------ */ /// \todo protected: does not compile with intel check why public: template ::value > class iterator_internal; public: /* ------------------------------------------------------------------------ */ // template class iterator : public iterator_internal {}; // template class const_iterator : public iterator_internal {}; /* ------------------------------------------------------------------------ */ //template using iterator = iterator_internal; template class iterator : public iterator_internal { public: typedef iterator_internal parent; typedef typename parent::value_type value_type; typedef typename parent::pointer pointer; typedef typename parent::reference reference; typedef typename parent::difference_type difference_type; typedef typename parent::iterator_category iterator_category; public: iterator() : parent() {}; iterator(pointer_type data, UInt offset) : parent(data, offset) {}; iterator(pointer warped) : parent(warped) {}; iterator(const iterator & it) : parent(it) {}; iterator(const parent & it) : parent(it) {}; inline iterator operator+(difference_type n) { return parent::operator+(n);; } inline iterator operator-(difference_type n) { return parent::operator-(n);; } inline difference_type operator-(const iterator & b) { return parent::operator-(b); } inline iterator & operator++() { parent::operator++(); return *this; }; inline iterator & operator--() { parent::operator--(); return *this; }; inline iterator & operator+=(const UInt n) { parent::operator+=(n); return *this; } }; /* ------------------------------------------------------------------------ */ //template using const_iterator = iterator_internal; template class const_iterator : public iterator_internal { public: typedef iterator_internal parent; typedef typename parent::value_type value_type; typedef typename parent::pointer pointer; typedef typename parent::reference reference; typedef typename parent::difference_type difference_type; typedef typename parent::iterator_category iterator_category; public: const_iterator() : parent() {}; const_iterator(pointer_type data, UInt offset) : parent(data, offset) {}; const_iterator(pointer warped) : parent(warped) {}; const_iterator(const parent & it) : parent(it) {}; const_iterator(const iterator & it) : parent(it) {}; inline const_iterator operator+(difference_type n) { return parent::operator+(n); } inline const_iterator operator-(difference_type n) { return parent::operator-(n); } inline difference_type operator-(const const_iterator & b) { return parent::operator-(b); } inline const_iterator & operator++() { parent::operator++(); return *this; }; inline const_iterator & operator--() { parent::operator--(); return *this; }; inline const_iterator & operator+=(const UInt n) { parent::operator+=(n); return *this; } }; inline iterator begin(); inline iterator end(); inline const_iterator begin() const; inline const_iterator end() const; inline iterator< Vector > begin(UInt n); inline iterator< Vector > end(UInt n); inline const_iterator< Vector > begin(UInt n) const; inline const_iterator< Vector > end(UInt n) const; inline iterator< Matrix > begin(UInt m, UInt n); inline iterator< Matrix > end(UInt m, UInt n); inline const_iterator< Matrix > begin(UInt m, UInt n) const; inline const_iterator< Matrix > end(UInt m, UInt n) const; /// /!\ to use with caution inline iterator< Vector > begin_reinterpret(UInt n, UInt size); inline iterator< Vector > end_reinterpret(UInt n, UInt size); inline const_iterator< Vector > begin_reinterpret(UInt n, UInt size) const; inline const_iterator< Vector > end_reinterpret(UInt n, UInt size) const; inline iterator< Matrix > begin_reinterpret(UInt m, UInt n, UInt size); inline iterator< Matrix > end_reinterpret(UInt m, UInt n, UInt size); inline const_iterator< Matrix > begin_reinterpret(UInt m, UInt n, UInt size) const; inline const_iterator< Matrix > end_reinterpret(UInt m, UInt n, UInt size) const; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// add an element at the end of the vector with the value value for all /// component inline void push_back(const_reference value); /// add an element at the end of the vector inline void push_back(const value_type new_elem[]); template inline void push_back(const iterator & it); inline void push_back(const Vector new_elem); /** * remove an element and move the last one in the hole * /!\ change the order in the vector */ inline void erase(UInt i); template inline void erase(const iterator & it); /// change the size of the vector and allocate more memory if needed void resize(UInt size); /// change the number of components by interlacing data void extendComponentsInterlaced(UInt multiplicator, UInt stride); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; // Array& operator=(const Array& vect); /// search elem in the vector, return the position of the first occurrence or /// -1 if not found Int find(const_reference elem) const; Int find(T elem[]) const; /// set a vvector to 0 inline void clear() { std::fill_n(values, size*nb_component, T()); }; + /// set a vector to the value t + inline void set(T & t) { std::fill_n(values, size*nb_component, t); }; + /// copy the content of an other vector void copy(const Array & vect); /// give the address of the memory allocated for this vector T * storage() const { return values; }; protected: /// perform the allocation for the constructors void allocate(UInt size, UInt nb_component = 1); /// resize without initializing the memory void resizeUnitialized(UInt new_size); /* ------------------------------------------------------------------------ */ /* Operators */ /* ------------------------------------------------------------------------ */ public: Array & operator-=(const Array & vect); Array & operator+=(const Array & vect); Array & operator*=(const T & alpha); inline reference operator()(UInt i, UInt j = 0); inline const_reference operator()(UInt i, UInt j = 0) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: UInt getSize() const{ return this->size; }; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ public: /// array of values T * values; // /!\ very dangerous }; __END_AKANTU__ #include "aka_types.hh" __BEGIN_AKANTU__ #include "aka_vector_tmpl.hh" /* -------------------------------------------------------------------------- */ /* Inline Functions Array */ /* -------------------------------------------------------------------------- */ template inline std::ostream & operator<<(std::ostream & stream, const Array & _this) { _this.printself(stream); return stream; } /* -------------------------------------------------------------------------- */ /* Inline Functions ArrayBase */ /* -------------------------------------------------------------------------- */ inline std::ostream & operator<<(std::ostream & stream, const ArrayBase & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_VECTOR_HH__ */ diff --git a/src/fem/element_class_tmpl.hh b/src/fem/element_class_tmpl.hh index af3b4b225..da370e1b7 100644 --- a/src/fem/element_class_tmpl.hh +++ b/src/fem/element_class_tmpl.hh @@ -1,397 +1,397 @@ /** * @file element_class_tmpl.hh * * @author Nicolas Richart * * @date Wed Nov 14 16:40:51 2012 * * @brief Implementation of the inline templated function of the element class * descriptions * * @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 . * */ /* -------------------------------------------------------------------------- */ /* GaussIntegrationElement */ /* -------------------------------------------------------------------------- */ template const Matrix GaussIntegrationElement::getQuadraturePoints() { Matrix quads(quad, - ElementClass::getNaturalSpaceDimension(), - nb_quadrature_points); + ElementClass::getNaturalSpaceDimension(), + nb_quadrature_points); return quads; } /* -------------------------------------------------------------------------- */ template const Vector GaussIntegrationElement::getWeights() { return Vector(weights, nb_quadrature_points); } /* -------------------------------------------------------------------------- */ /* GeometricalElement */ /* -------------------------------------------------------------------------- */ template inline const Matrix GeometricalElement::getFacetLocalConnectivityPerElement() { return Matrix(facet_connectivity, nb_facets, nb_nodes_per_facet); } /* -------------------------------------------------------------------------- */ template inline bool GeometricalElement::contains(const Vector & coords) { return GeometricalShapeContains::contains(coords); } /* -------------------------------------------------------------------------- */ template<> inline bool GeometricalShapeContains<_gst_point>::contains(const Vector & coords) { return (coords(0) < std::numeric_limits::epsilon()); } /* -------------------------------------------------------------------------- */ template<> inline bool GeometricalShapeContains<_gst_square>::contains(const Vector & coords) { bool in = true; for (UInt i = 0; i < coords.size() && in; ++i) in &= ((coords(i) >= -1.) && (coords(i) <= 1.)); return in; } /* -------------------------------------------------------------------------- */ template<> inline bool GeometricalShapeContains<_gst_triangle>::contains(const Vector & coords) { bool in = true; Real sum = 0; for (UInt i = 0; (i < coords.size()) && in; ++i) { in &= ((coords(i) >= 0) && (coords(i) <= 1.)); sum += coords(i); } if(in) return (in && (sum <= 1)); return in; } /* -------------------------------------------------------------------------- */ /* 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); computeShapes(natural_coord(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); computeDNDS(natural_coord(p), dnds_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(nb_nodes_per_element); computeShapes(natural_coords, shapes); Matrix interpm(interpolated.storage(), nodal_values.rows(), 1); Matrix shapesm(shapes.storage(), nb_nodes_per_element, 1); interpm.mul(nodal_values, shapesm); } /* -------------------------------------------------------------------------- */ /// @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(natural_space_dimension, 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); 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_element::natural_space_dimension, interpolation_element::nb_nodes_per_element); Matrix J(natural_coords.rows(), node_coords.rows()); for (UInt p = 0; p < nb_points; ++p) { interpolation_element::computeDNDS(natural_coords(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_element::natural_space_dimension, "cannot extract a normal because of dimension mismatch " << dimension << " " << interpolation_element::natural_space_dimension); Matrix J(dimension, interpolation_element::natural_space_dimension); for (UInt p = 0; p < nb_points; ++p) { interpolation_element::gradientOnNaturalCoordinates(coord(p), f, J); if (dimension == 2) { Math::normal2(J.storage(), normals(p).storage()); } if (dimension == 3){ Math::normal3(J(0).storage(), J(1).storage(), normals(p).storage()); } } } /* ------------------------------------------------------------------------- */ /** * In the non linear cases we need to iterate to find the natural coordinates @f$\xi@f$ * provided real coordinates @f$x@f$. * * We want to solve: @f$ x- \phi(\xi) = 0@f$ with @f$\phi(\xi) = \sum_I N_I(\xi) x_I@f$ * the mapping function which uses the nodal coordinates @f$x_I@f$. * * To that end we use the Newton method and the following series: * * @f$ \frac{\partial \phi(x_k)}{\partial \xi} \left( \xi_{k+1} - \xi_k \right) = x - \phi(x_k)@f$ * * When we consider elements embedded in a dimension higher than them (2D triangle in a 3D space for example) * @f$ J = \frac{\partial \phi(\xi_k)}{\partial \xi}@f$ is of dimension @f$dim_{space} \times dim_{elem}@f$ which * is not invertible in most cases. Rather we can solve the problem: * * @f$ J^T J \left( \xi_{k+1} - \xi_k \right) = J^T \left( x - \phi(\xi_k) \right) @f$ * * So that * * @f$ d\xi = \xi_{k+1} - \xi_k = (J^T J)^{-1} J^T \left( x - \phi(\xi_k) \right) @f$ * * So that if the series converges we have: * * @f$ 0 = J^T \left( \phi(\xi_\infty) - x \right) @f$ * * And we see that this is ill-posed only if @f$ J^T x = 0@f$ which means that the vector provided * is normal to any tangent which means it is outside of the element itself. * * @param real_coords: the real coordinates the natural coordinates are sought for * @param node_coords: the coordinates of the nodes forming the element * @param natural_coords: output->the sought natural coordinates * @param spatial_dimension: spatial dimension of the problem * **/ template inline void ElementClass::inverseMap(const Vector & real_coords, const Matrix & node_coords, Vector & natural_coords, Real tolerance) { UInt spatial_dimension = real_coords.size(); UInt dimension = natural_coords.size(); //matrix copy of the real_coords Matrix mreal_coords(real_coords.storage(), spatial_dimension, 1); //initial guess // Matrix natural_guess(natural_coords.storage(), dimension, 1); natural_coords.clear(); // real space coordinates provided by initial guess Matrix physical_guess(dimension, 1); // objective function f = real_coords - physical_guess Matrix f(dimension, 1); // dnds computed on the natural_guess // Matrix dnds(interpolation_element::nb_nodes_per_element, spatial_dimension); // J Jacobian matrix computed on the natural_guess Matrix J(spatial_dimension, dimension); // G = J^t * J Matrix G(spatial_dimension, spatial_dimension); // Ginv = G^{-1} Matrix Ginv(spatial_dimension, spatial_dimension); // J = Ginv * J^t Matrix F(spatial_dimension, dimension); // dxi = \xi_{k+1} - \xi in the iterative process Matrix dxi(spatial_dimension, 1); /* --------------------------- */ /* init before iteration loop */ /* --------------------------- */ // do interpolation Vector physical_guess_v(physical_guess.storage(), dimension); interpolation_element::interpolateOnNaturalCoordinates(natural_coords, node_coords, physical_guess_v); // compute initial objective function value f = real_coords - physical_guess f = mreal_coords; f -= physical_guess; // compute initial error Real inverse_map_error = f.norm(); /* --------------------------- */ /* iteration loop */ /* --------------------------- */ while(tolerance < inverse_map_error) { //compute J^t interpolation_element::gradientOnNaturalCoordinates(natural_coords, node_coords, J); //compute G G.mul(J, J); // inverse G Ginv.inverse(G); //compute F F.mul(Ginv, J); //compute increment dxi.mul(F, f); //update our guess natural_coords += dxi(0); //interpolate interpolation_element::interpolateOnNaturalCoordinates(natural_coords, node_coords, physical_guess_v); // compute error f = mreal_coords; f -= physical_guess; inverse_map_error = f.norm(); } // memcpy(natural_coords.storage(), natural_guess.storage(), sizeof(Real) * natural_coords.size()); } diff --git a/src/fem/fem_template_tmpl.hh b/src/fem/fem_template_tmpl.hh index 6d96f8e6a..5ba9b5783 100644 --- a/src/fem/fem_template_tmpl.hh +++ b/src/fem/fem_template_tmpl.hh @@ -1,924 +1,919 @@ /** * @file fem_template_tmpl.hh * @author Nicolas Richart * @date Mon Nov 5 17:08:50 2012 * * @brief Template implementation of FEMTemplate * * @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 . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template