diff --git a/CMakeLists.txt b/CMakeLists.txt index e527c0c84..222803a07 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,341 +1,351 @@ #=============================================================================== # @file CMakeLists.txt # @author Nicolas Richart # @date Fri Jun 11 09:46:59 2010 # # @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 # #=============================================================================== #=============================================================================== # _ ,-, _ # ,--, /: :\/': :`\/: :\ # |`; ' `,' `.; `: | _ _ # | | | ' | |. | | | | # | : | F E | A S | T++ || __ _| | ____ _ _ __ | |_ _ _ # | :. | : | : | : | \ / _` | |/ / _` | '_ \| __| | | | # \__/: :.. : :.. | :.. | ) | (_| | < (_| | | | | |_| |_| | # `---',\___/,\___/ /' \__,_|_|\_\__,_|_| |_|\__|\__,_| # `==._ .. . /' # `-::-' #=============================================================================== #=============================================================================== # CMake Project #=============================================================================== cmake_minimum_required(VERSION 2.6) project(AKANTU) enable_language(CXX) #=============================================================================== # Misc. #=============================================================================== set(AKANTU_CMAKE_DIR "${AKANTU_SOURCE_DIR}/cmake") set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake") set(BUILD_SHARED_LIBS ON CACHE BOOL "Build shared libraries.") #=============================================================================== # Version Number #=============================================================================== # AKANTU version number. An even minor number corresponds to releases. set(AKANTU_MAJOR_VERSION 0) set(AKANTU_MINOR_VERSION 1) set(AKANTU_BUILD_VERSION 0) set(AKANTU_VERSION "${AKANTU_MAJOR_VERSION}.${AKANTU_MINOR_VERSION}.${AKANTU_BUILD_VERSION}" ) # Append the library version information to the library target properties if(NOT AKANTU_NO_LIBRARY_VERSION) set(AKANTU_LIBRARY_PROPERTIES ${AKANTU_LIBRARY_PROPERTIES} VERSION "${AKANTU_VERSION}" SOVERSION "${AKANTU_MAJOR_VERSION}.${AKANTU_MINOR_VERSION}" ) endif(NOT AKANTU_NO_LIBRARY_VERSION) #=============================================================================== # Options #=============================================================================== option(AKANTU_DEBUG "Compiles akantu with the debug messages" ON) option(AKANTU_DOCUMENTATION "Build source documentation using Doxygen." OFF) # compilation switch option(AKANTU_BUILD_CONTACT "Build the contact algorithm" OFF) # features option(AKANTU_USE_IOHELPER "Add IOHelper support in akantu" OFF) option(AKANTU_USE_MPI "Add MPI support in akantu" OFF) option(AKANTU_USE_SCOTCH "Add Scotch support in akantu" OFF) option(AKANTU_USE_MUMPS "Add Mumps support in akantu" OFF) option(AKANTU_USE_BLAS "Use BLAS for arithmetic operations" OFF) option(AKANTU_USE_EPSN "Use EPSN streering environment" OFF) #=============================================================================== # Options handling #=============================================================================== +#=========================================================================== +# Documentation +#=========================================================================== +find_package(Boost REQUIRED) +if(BOOST_FOUND) + include_directories(BOOST_INCLUDE_PATH) +endif(BOOST_FOUND) + + # Debug if(NOT AKANTU_DEBUG) add_definitions(-DAKANTU_NDEBUG) endif(NOT AKANTU_DEBUG) # IOHelper if(AKANTU_USE_IOHELPER) find_package(IOHelper REQUIRED) if(IOHELPER_FOUND) add_definitions(-DAKANTU_USE_IOHELPER) set(AKANTU_EXTERNAL_LIB_INCLUDE_PATH ${AKANTU_EXTERNAL_LIB_INCLUDE_PATH} ${IOHELPER_INCLUDE_PATH} ) set(AKANTU_EXTERNAL_LIBRARIES ${AKANTU_EXTERNAL_LIBRARIES} ${IOHELPER_LIBRARIES} ) endif(IOHELPER_FOUND) endif(AKANTU_USE_IOHELPER) # MPI set(AKANTU_MPI_ON FALSE) if(AKANTU_USE_MPI) find_package(MPI REQUIRED) if(MPI_FOUND) add_definitions(-DAKANTU_USE_MPI) set(AKANTU_EXTERNAL_LIB_INCLUDE_PATH ${AKANTU_EXTERNAL_LIB_INCLUDE_PATH} ${MPI_INCLUDE_PATH} ) set(AKANTU_EXTERNAL_LIBRARIES ${AKANTU_EXTERNAL_LIBRARIES} ${MPI_LIBRARY} ${MPI_EXTRA_LIBRARY} ) set(AKANTU_MPI_ON TRUE) endif(MPI_FOUND) endif(AKANTU_USE_MPI) # Scotch set(AKANTU_SCOTCH_ON FALSE) if(AKANTU_USE_SCOTCH) find_package(Scotch REQUIRED) if(SCOTCH_FOUND) add_definitions(-DAKANTU_USE_SCOTCH) set(AKANTU_EXTERNAL_LIB_INCLUDE_PATH ${AKANTU_EXTERNAL_LIB_INCLUDE_PATH} ${SCOTCH_INCLUDE_PATH} ) set(AKANTU_EXTERNAL_LIBRARIES ${AKANTU_EXTERNAL_LIBRARIES} ${SCOTCH_LIBRARIES} ) set(AKANTU_SCOTCH_ON TRUE) endif(SCOTCH_FOUND) endif(AKANTU_USE_SCOTCH) # MUMPS set(AKANTU_MUMPS_ON FALSE) if(AKANTU_USE_MUMPS) find_package(Mumps REQUIRED) if(MUMPS_FOUND) add_definitions(-DAKANTU_USE_MUMPS) set(AKANTU_EXTERNAL_LIB_INCLUDE_PATH ${AKANTU_EXTERNAL_LIB_INCLUDE_PATH} ${MUMPS_INCLUDE_PATH} ) set(AKANTU_EXTERNAL_LIBRARIES ${AKANTU_EXTERNAL_LIBRARIES} ${MUMPS_LIBRARIES} ) set(AKANTU_MUMPS_ON TRUE) endif(MUMPS_FOUND) endif(AKANTU_USE_MUMPS) # BLAS set(AKANTU_BLAS_ON FALSE) if(AKANTU_USE_BLAS) enable_language(Fortran) find_package(BLAS) if(BLAS_FOUND) add_definitions(-DAKANTU_USE_CBLAS) set(AKANTU_EXTERNAL_LIBRARIES ${AKANTU_EXTERNAL_LIBRARIES} ${BLAS_LIBRARIES} ) set(AKANTU_BLAS_ON TRUE) endif(BLAS_FOUND) endif(AKANTU_USE_BLAS) # EPSN set(AKANTU_EPSN_ON FALSE) if(AKANTU_USE_EPSN) find_package(EPSN) if(EPSN_FOUND) add_definitions(-DAKANTU_USE_EPSN) set(AKANTU_EXTERNAL_LIB_INCLUDE_PATH ${AKANTU_EXTERNAL_LIB_INCLUDE_PATH} ${EPSN_INCLUDE_DIR} ) set(AKANTU_EXTERNAL_LIBRARIES ${AKANTU_EXTERNAL_LIBRARIES} ${EPSN_LIBRARIES} ) set(AKANTU_EPSN_ON TRUE) endif(EPSN_FOUND) endif(AKANTU_USE_EPSN) #=============================================================================== # Library #=============================================================================== set(AKANTU_COMMON_SRC common/aka_common.cc common/aka_error.cc common/aka_extern.cc common/aka_static_memory.cc common/aka_memory.cc common/aka_vector.cc common/aka_math.cc fem/mesh.cc fem/fem.cc fem/element_class.cc # model/integration_scheme/central_difference.cc model/solid_mechanics/solid_mechanics_model.cc model/solid_mechanics/solid_mechanics_model_mass.cc model/solid_mechanics/solid_mechanics_model_boundary.cc model/solid_mechanics/solid_mechanics_model_material.cc model/solid_mechanics/material.cc model/solid_mechanics/materials/material_elastic.cc mesh_utils/mesh_io.cc + mesh_utils/mesh_pbc.cc mesh_utils/mesh_io/mesh_io_msh.cc mesh_utils/mesh_partition.cc mesh_utils/mesh_utils.cc solver/sparse_matrix.cc solver/solver.cc synchronizer/ghost_synchronizer.cc synchronizer/synchronizer.cc synchronizer/communicator.cc synchronizer/static_communicator.cc ) set(AKANTU_CONTACT_SRC model/solid_mechanics/contact.cc model/solid_mechanics/contact_search.cc model/solid_mechanics/contact_neighbor_structure.cc model/solid_mechanics/contact/contact_2d_explicit.cc model/solid_mechanics/contact/contact_search_2d_explicit.cc model/solid_mechanics/contact/regular_grid_neighbor_structure.cc model/solid_mechanics/contact/contact_search_3d_explicit.cc model/solid_mechanics/contact/contact_3d_explicit.cc model/solid_mechanics/contact/grid_2d_neighbor_structure.cc model/solid_mechanics/contact/contact_rigid_no_friction.cc ) if(AKANTU_BUILD_CONTACT) set(AKANTU_COMMON_SRC ${AKANTU_COMMON_SRC} ${AKANTU_CONTACT_SRC}) endif(AKANTU_BUILD_CONTACT) if(AKANTU_USE_MPI AND MPI_FOUND) set(AKANTU_COMMON_SRC ${AKANTU_COMMON_SRC} synchronizer/static_communicator_mpi.cc ) endif(AKANTU_USE_MPI AND MPI_FOUND) if(AKANTU_SCOTCH_ON) set(AKANTU_COMMON_SRC ${AKANTU_COMMON_SRC} mesh_utils/mesh_partition/mesh_partition_scotch.cc ) endif(AKANTU_SCOTCH_ON) if(AKANTU_MUMPS_ON) set(AKANTU_COMMON_SRC ${AKANTU_COMMON_SRC} solver/solver_mumps.cc ) endif(AKANTU_MUMPS_ON) set(AKANTU_INCLUDE_DIRS common fem/ mesh_utils/ mesh_utils/mesh_io/ mesh_utils/mesh_partition/ model/ model/integration_scheme model/solid_mechanics model/solid_mechanics/materials model/solid_mechanics/contact synchronizer/ solver/ ) include_directories( ${AKANTU_INCLUDE_DIRS} ${AKANTU_EXTERNAL_LIB_INCLUDE_PATH} ) add_library(akantu ${AKANTU_COMMON_SRC}) set_target_properties(akantu PROPERTIES ${AKANTU_LIBRARY_PROPERTIES}) #=========================================================================== # Documentation #=========================================================================== if(AKANTU_DOCUMENTATION) find_package(Doxygen REQUIRED) if(DOXYGEN_FOUND) set(DOXYGEN_WARNINGS NO) set(DOXYGEN_QUIET YES) if(CMAKE_VERBOSE_MAKEFILE) set(DOXYGEN_WARNINGS YES) set(DOXYGEN_QUIET NO) endif(CMAKE_VERBOSE_MAKEFILE) add_subdirectory(doc/doxygen) endif(DOXYGEN_FOUND) endif(AKANTU_DOCUMENTATION) #=============================================================================== # Tests #=============================================================================== ENABLE_TESTING() INCLUDE(CTest) option(AKANTU_TESTS "Activate tests" OFF) if(AKANTU_TESTS) find_package(GMSH REQUIRED) add_subdirectory(test) endif(AKANTU_TESTS) #=============================================================================== # Examples #=============================================================================== option(AKANTU_EXAMPLES "Activate examples" OFF) if(AKANTU_EXAMPLES) find_package(GMSH REQUIRED) add_subdirectory(examples) endif(AKANTU_EXAMPLES) diff --git a/common/aka_common.hh b/common/aka_common.hh index 64f7cdc04..9ecd817a9 100644 --- a/common/aka_common.hh +++ b/common/aka_common.hh @@ -1,295 +1,326 @@ /** * @file aka_common.hh * @author Nicolas Richart * @date Fri Jun 11 09:48:06 2010 * * @namespace 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 #include /* -------------------------------------------------------------------------- */ #include #include #include #include #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ #define __BEGIN_AKANTU__ namespace akantu { #define __END_AKANTU__ } /* -------------------------------------------------------------------------- */ #include "aka_error.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* Common types */ /* -------------------------------------------------------------------------- */ typedef double Real; typedef unsigned int UInt; typedef int Int; typedef std::string ID; #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; typedef ID VectorID; /* -------------------------------------------------------------------------- */ /* Mesh/FEM/Model types */ /* -------------------------------------------------------------------------- */ typedef ID MeshID; typedef ID FEMID; typedef ID ModelID; typedef ID MaterialID; typedef ID SparseMatrixID; typedef ID SolverID; typedef UInt Surface; + +/* -------------------------------------------------------------------------- */ +// 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(macro) \ + switch(type) { \ + BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_CASE_MACRO,macro,AKANTU_ELEMENT_TYPE) \ + case _not_defined: \ + case _max_element_type: { \ + AKANTU_DEBUG_ERROR("Wrong type : " << type); \ + break; \ + } \ + } +/* -------------------------------------------------------------------------- */ + +/// @boost sequence of element to loop on in global tasks +#define AKANTU_ELEMENT_TYPE \ + (_segment_2) \ + (_segment_3) \ + (_triangle_3) \ + (_triangle_6) \ + (_tetrahedron_4) \ + (_tetrahedron_10) \ + (_quadrangle_4) \ + (_point) + + /// @enum ElementType type of element potentially contained in a Mesh enum ElementType { _not_defined = 0, _segment_2 = 1, /// first order segment _segment_3 = 2, /// second order segment _triangle_3 = 3, /// first order triangle _triangle_6 = 4, /// second order triangle _tetrahedron_4 = 5, /// first order tetrahedron _tetrahedron_10 = 6, /// second order tetrahedron @remark not implemented yet _quadrangle_4, /// first order quadrangle _max_element_type, _point /// point only for some algorithm to be generic like mesh partitioning }; /// @enum MaterialType different materials implemented enum MaterialType { _elastic = 0, _max_material_type }; typedef void (*BoundaryFunction)(double *,double *); /// @enum BoundaryFunctionType type of function passed for boundary conditions enum BoundaryFunctionType { _bft_stress, _bft_forces }; /// @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_no_fric = 3 }; enum ContactSearchType { _cst_not_defined = 0, _cst_2d_expli = 1, _cst_3d_expli = 2 }; enum ContactNeighborStructureType { _cnst_not_defined = 0, _cnst_regular_grid = 1, _cnst_2d_grid = 2 }; /* -------------------------------------------------------------------------- */ /* Ghosts handling */ /* -------------------------------------------------------------------------- */ typedef ID SynchronizerID; /// @enum CommunicatorType type of communication method to use enum CommunicatorType { _communicator_mpi, _communicator_dummy }; /// @enum GhostSynchronizationTag type of synchronizations enum GhostSynchronizationTag { /// 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 /// Test tag _gst_test }; /// @enum GhostType type of ghost enum GhostType { _not_ghost, _ghost, _casper // not used but a real cute ghost }; /// @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_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_ELEMENT_TYPE(name, variable, type) \ inline type get##name (const ::akantu::ElementType & el_type) const { \ AKANTU_DEBUG_IN(); \ AKANTU_DEBUG_ASSERT(variable[el_type] != NULL, \ "No " << #variable << " of type " \ << el_type << " in " << id); \ AKANTU_DEBUG_OUT(); \ return *variable[el_type]; \ } /* -------------------------------------------------------------------------- */ //! standard output stream operator for ElementType inline std::ostream & operator <<(std::ostream & stream, ElementType type) { switch(type) { case _segment_2 : stream << "segment_2" ; break; case _segment_3 : stream << "segment_3" ; break; case _triangle_3 : stream << "triangle_3" ; break; case _triangle_6 : stream << "triangle_6" ; break; case _tetrahedron_4 : stream << "tetrahedron_4" ; break; case _tetrahedron_10 : stream << "tetrahedron_10"; break; case _quadrangle_4 : stream << "quadrangle_4" ; break; case _not_defined : stream << "undefined" ; break; case _max_element_type : stream << "ElementType(" << (int) type << ")"; break; case _point : stream << "point"; break; } return stream; } /* -------------------------------------------------------------------------- */ 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 void to_lower(std::string & str) { std::transform(str.begin(), str.end(), str.begin(), (int(*)(int))std::tolower); } /* -------------------------------------------------------------------------- */ inline void trim(std::string & to_trim) { size_t first = to_trim.find_first_not_of(" \t"); if (first != std::string::npos) { size_t last = to_trim.find_last_not_of(" \t"); to_trim = to_trim.substr(first, last - first + 1); } else to_trim = ""; } __END_AKANTU__ #endif /* __AKANTU_COMMON_HH__ */ diff --git a/common/aka_extern.cc b/common/aka_extern.cc index fe40834ad..0a23ecb36 100644 --- a/common/aka_extern.cc +++ b/common/aka_extern.cc @@ -1,62 +1,66 @@ /** * @file aka_extern.cc * @author Nicolas Richart * @date Mon Jun 14 14:34:14 2010 * * @brief initialisation of all global variables * to insure the order of creation * * @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 /* -------------------------------------------------------------------------- */ #include "aka_common.hh" +#include "aka_math.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /** \todo write function to get this * values from the environment or a config file */ /* -------------------------------------------------------------------------- */ /* error.hpp variables */ /* -------------------------------------------------------------------------- */ namespace debug { /// standard output for debug messages std::ostream & _akantu_debug_cout = std::cerr; /// standard output for normal messages std::ostream & _akantu_cout = std::cout; /// debug level DebugLevel _debug_level = (DebugLevel) 5; /// parallel context used in debug messages std::string _parallel_context = ""; + } +Real Math::tolerance = 1e-8; + /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/common/aka_math.hh b/common/aka_math.hh index 22d4601b6..2dd523f28 100644 --- a/common/aka_math.hh +++ b/common/aka_math.hh @@ -1,183 +1,188 @@ /** * @file aka_math.hh * @author Nicolas Richart * @date Wed Jul 28 11:51:56 2010 * * @brief mathematical operations * * @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_MATH_H__ #define __AKANTU_AKA_MATH_H__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_vector.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ class Math { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Matrix algebra */ /* ------------------------------------------------------------------------ */ /// @f$ y = A*x @f$ static void matrix_vector(UInt m, UInt n, const Vector & A, const Vector & x, Vector & y); /// @f$ y = A*x @f$ static inline void matrix_vector(UInt m, UInt n, const Real * A, const Real * x, Real * y); /// @f$ C = A*B @f$ static void matrix_matrix(UInt m, UInt n, UInt k, const Vector & A, const Vector & B, Vector & C); /// @f$ C = A*B^t @f$ static void matrix_matrixt(UInt m, UInt n, UInt k, const Vector & A, const Vector & B, Vector & C); /// @f$ C = A*B @f$ static inline void matrix_matrix(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C); /// @f$ C = A^t*B @f$ static inline void matrixt_matrix(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C); /// @f$ C = A*B^t @f$ static inline void matrix_matrixt(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C); /// @f$ C = A^t*B^t @f$ static inline void matrixt_matrixt(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C); /// determinent of a 3x3 matrix static inline Real det3(const Real * mat); /// determinent of a 2x2 matrix static inline Real det2(const Real * mat); /// inverse a 3x3 matrix static inline void inv3(const Real * mat, Real * inv); /// inverse a 2x2 matrix static inline void inv2(const Real * mat, Real * inv); /// vector cross product static inline void vectorProduct3(const Real * v1, const Real * v2, Real * res); /// compute normal a normal to a vector static inline void normal2(const Real * v1, Real * res); /// compute normal a normal to a vector static inline void normal3(const Real * v1,const Real * v2, Real * res); /// normalize a vector static inline void normalize2(Real * v); /// normalize a vector static inline void normalize3(Real * v); /// return norm of a 2-vector static inline Real norm2(const Real * v); /// return norm of a 3-vector static inline Real norm3(const Real * v); /// return the dot product between 2 vectors in 2d static inline Real vectorDot2(const Real * v1, const Real * v2); /// return the dot product between 2 vectors in 3d static inline Real vectorDot3(const Real * v1, const Real * v2); /* ------------------------------------------------------------------------ */ /* Geometry */ /* ------------------------------------------------------------------------ */ /// distance in 2D between x and y static inline Real distance_2d(const Real * x, const Real * y); /// distance in 3D between x and y static inline Real distance_3d(const Real * x, const Real * y); /// radius of the in-circle of a triangle static inline Real triangle_inradius(const Real * coord1, const Real * coord2, const Real * coord3); /// radius of the in-circle of a tetrahedron static inline Real tetrahedron_inradius(const Real * coord1, const Real * coord2, const Real * coord3, const Real * coord4); /// volume of a tetrahedron static inline Real tetrahedron_volume(const Real * coord1, const Real * coord2, const Real * coord3, const Real * coord4); /// compute the barycenter of n points static inline void barycenter(const Real * coord, UInt nb_points, UInt spatial_dimension, Real * barycenter); /// vector between x and y static inline void vector_2d(const Real * x, const Real * y, Real * vec); /// vector pointing from x to y in 3 spatial dimension static inline void vector_3d(const Real * x, const Real * y, Real * vec); + /// test if two scalar are equal within a given tolerance + static inline bool are_float_equal(const Real x, const Real y); + + /// tolerance for functions that need one + static Real tolerance; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "aka_math_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_AKA_MATH_H__ */ diff --git a/common/aka_math_inline_impl.cc b/common/aka_math_inline_impl.cc index 02154fa8e..b3cf1dc7b 100644 --- a/common/aka_math_inline_impl.cc +++ b/common/aka_math_inline_impl.cc @@ -1,387 +1,391 @@ /** * @file aka_math_inline_impl.cc * @author Nicolas Richart * @date Wed Jul 28 13:20:35 2010 * * @brief Implementation of the inline functions of the math toolkit * * @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 . * */ #ifdef AKANTU_USE_CBLAS # ifndef AKANTU_USE_CBLAS_MKL # include # else // AKANTU_USE_CBLAS_MKL # include # endif //AKANTU_USE_CBLAS_MKL #endif //AKANTU_USE_CBLAS /* -------------------------------------------------------------------------- */ inline void Math::matrix_vector(UInt m, UInt n, const Real * A, const Real * x, Real * y) { #ifdef AKANTU_USE_CBLAS /// y = alpha*op(A)*x + beta*y cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n, 1, A, n, x, 1, 0, y, 1); #else memset(y, 0, m*sizeof(Real)); for (UInt i = 0; i < m; ++i) { UInt A_i = i * n; for (UInt j = 0; j < n; ++j) { y[i] += A[A_i + j] * x[j]; } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrix_matrix(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C) { #ifdef AKANTU_USE_CBLAS /// C := alpha*op(A)*op(B) + beta*C cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, 1, A, k, B, n, 0, C, n); #else memset(C, 0, m*n*sizeof(Real)); for (UInt i = 0; i < m; ++i) { UInt A_i = i * k; UInt C_i = i * n; for (UInt j = 0; j < n; ++j) { for (UInt l = 0; l < k; ++l) { C[C_i + j] += A[A_i + l] * B[l * n + j]; } } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrixt_matrix(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C) { #ifdef AKANTU_USE_CBLAS /// C := alpha*op(A)*op(B) + beta*C cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans, m, n, k, 1, A, m, B, n, 0, C, n); #else memset(C, 0, m*n*sizeof(Real)); for (UInt i = 0; i < m; ++i) { // UInt A_i = i * k; UInt C_i = i * n; for (UInt j = 0; j < n; ++j) { for (UInt l = 0; l < k; ++l) { C[C_i + j] += A[l * m + i] * B[l * n + j]; } } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrix_matrixt(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C) { #ifdef AKANTU_USE_CBLAS /// C := alpha*op(A)*op(B) + beta*C cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, m, n, k, 1, A, k, B, k, 0, C, n); #else memset(C, 0, m*n*sizeof(Real)); for (UInt i = 0; i < m; ++i) { UInt A_i = i * k; UInt C_i = i * n; for (UInt j = 0; j < n; ++j) { UInt B_j = j * k; for (UInt l = 0; l < k; ++l) { C[C_i + j] += A[A_i + l] * B[B_j + l]; } } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrixt_matrixt(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C) { #ifdef AKANTU_USE_CBLAS /// C := alpha*op(A)*op(B) + beta*C cblas_dgemm(CblasRowMajor, CblasTrans, CblasTrans, m, n, k, 1, A, m, B, k, 0, C, n); #else memset(C, 0, m * n * sizeof(Real)); for (UInt i = 0; i < m; ++i) { UInt C_i = i * n; for (UInt j = 0; j < n; ++j) { UInt B_j = j * k; for (UInt l = 0; l < k; ++l) { C[C_i + j] += A[l * m + i] * B[B_j + l]; } } } #endif } /* -------------------------------------------------------------------------- */ inline Real Math::det2(const Real * mat) { return mat[0]*mat[3] - mat[1]*mat[2]; } /* -------------------------------------------------------------------------- */ inline Real Math::det3(const Real * mat) { return mat[0]*(mat[4]*mat[8]-mat[7]*mat[5]) - mat[3]*(mat[1]*mat[8]-mat[7]*mat[2]) + mat[6]*(mat[1]*mat[5]-mat[4]*mat[2]); } /* -------------------------------------------------------------------------- */ inline void Math::normal2(const Real * vec,Real * normal) { normal[0] = vec[1]; normal[1] = -vec[0]; Math::normalize2(normal); } /* -------------------------------------------------------------------------- */ inline void Math::normal3(const Real * vec1,const Real * vec2,Real * normal) { Math::vectorProduct3(vec1,vec2,normal); Math::normalize3(normal); } /* -------------------------------------------------------------------------- */ inline void Math::normalize2(Real * vec) { Real norm = Math::norm2(vec); vec[0] /= norm; vec[1] /= norm; } /* -------------------------------------------------------------------------- */ inline void Math::normalize3(Real * vec) { Real norm = Math::norm3(vec); vec[0] /= norm; vec[1] /= norm; vec[2] /= norm; } /* -------------------------------------------------------------------------- */ inline Real Math::norm2(const Real * vec) { return sqrt(vec[0]*vec[0] + vec[1]*vec[1]); } /* -------------------------------------------------------------------------- */ inline Real Math::norm3(const Real * vec) { return sqrt(vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2]); } /* -------------------------------------------------------------------------- */ inline void Math::inv2(const Real * mat,Real * inv) { Real det_mat = det2(mat); inv[0] = mat[3] / det_mat; inv[1] = -mat[1] / det_mat; inv[2] = -mat[2] / det_mat; inv[3] = mat[0] / det_mat; } /* -------------------------------------------------------------------------- */ inline void Math::inv3(const Real * mat,Real * inv) { Real det_mat = det3(mat); inv[0] = (mat[4]*mat[8] - mat[7]*mat[5])/det_mat; inv[1] = (mat[2]*mat[7] - mat[8]*mat[1])/det_mat; inv[2] = (mat[1]*mat[5] - mat[4]*mat[2])/det_mat; inv[3] = (mat[5]*mat[6] - mat[8]*mat[3])/det_mat; inv[4] = (mat[0]*mat[8] - mat[6]*mat[2])/det_mat; inv[5] = (mat[2]*mat[3] - mat[5]*mat[0])/det_mat; inv[6] = (mat[3]*mat[7] - mat[6]*mat[4])/det_mat; inv[7] = (mat[1]*mat[6] - mat[7]*mat[0])/det_mat; inv[8] = (mat[0]*mat[4] - mat[3]*mat[1])/det_mat; } /* -------------------------------------------------------------------------- */ inline void Math::vectorProduct3(const Real * v1, const Real * v2, Real * res) { res[0] = v1[1]*v2[2] - v1[2]*v2[1]; res[1] = v1[2]*v2[0] - v1[0]*v2[2]; res[2] = v1[0]*v2[1] - v1[1]*v2[0]; } /* -------------------------------------------------------------------------- */ inline Real Math::vectorDot2(const Real * v1, const Real * v2) { return (v1[0]*v2[0] + v1[1]*v2[1]); } /* -------------------------------------------------------------------------- */ inline Real Math::vectorDot3(const Real * v1, const Real * v2) { return (v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]); } /* -------------------------------------------------------------------------- */ inline Real Math::distance_2d(const Real * x, const Real * y) { return sqrt((y[0] - x[0])*(y[0] - x[0]) + (y[1] - x[1])*(y[1] - x[1])); } /* -------------------------------------------------------------------------- */ inline Real Math::triangle_inradius(const Real * coord1, const Real * coord2, const Real * coord3) { /** * @f{eqnarray*}{ * r &=& A / s \\ * A &=& 1/4 * \sqrt{(a + b + c) * (a - b + c) * (a + b - c) (-a + b + c)} \\ * s &=& \frac{a + b + c}{2} * @f} */ Real a, b, c; a = distance_2d(coord1, coord2); b = distance_2d(coord2, coord3); c = distance_2d(coord1, coord3); Real s; s = (a + b + c) * 0.5; return sqrt((s - a) * (s - b) * (s - c) / s); } /* -------------------------------------------------------------------------- */ inline Real Math::distance_3d(const Real * x, const Real * y) { return sqrt((y[0] - x[0])*(y[0] - x[0]) + (y[1] - x[1])*(y[1] - x[1]) + (y[2] - x[2])*(y[2] - x[2]) ); } /* -------------------------------------------------------------------------- */ inline Real Math::tetrahedron_volume(const Real * coord1, const Real * coord2, const Real * coord3, const Real * coord4) { Real xx[9], vol; xx[0] = coord2[0]; xx[1] = coord2[1]; xx[2] = coord2[2]; xx[3] = coord3[0]; xx[4] = coord3[1]; xx[5] = coord3[2]; xx[6] = coord4[0]; xx[7] = coord4[1]; xx[8] = coord4[2]; vol = det3(xx); xx[0] = coord1[0]; xx[1] = coord1[1]; xx[2] = coord1[2]; xx[3] = coord3[0]; xx[4] = coord3[1]; xx[5] = coord3[2]; xx[6] = coord4[0]; xx[7] = coord4[1]; xx[8] = coord4[2]; vol -= det3(xx); xx[0] = coord1[0]; xx[1] = coord1[1]; xx[2] = coord1[2]; xx[3] = coord2[0]; xx[4] = coord2[1]; xx[5] = coord2[2]; xx[6] = coord4[0]; xx[7] = coord4[1]; xx[8] = coord4[2]; vol += det3(xx); xx[0] = coord1[0]; xx[1] = coord1[1]; xx[2] = coord1[2]; xx[3] = coord2[0]; xx[4] = coord2[1]; xx[5] = coord2[2]; xx[6] = coord3[0]; xx[7] = coord3[1]; xx[8] = coord3[2]; vol -= det3(xx); vol /= 6; return vol; } /* -------------------------------------------------------------------------- */ inline Real Math::tetrahedron_inradius(const Real * coord1, const Real * coord2, const Real * coord3, const Real * coord4) { Real l12, l13, l14, l23, l24, l34; l12 = distance_3d(coord1, coord2); l13 = distance_3d(coord1, coord3); l14 = distance_3d(coord1, coord4); l23 = distance_3d(coord2, coord3); l24 = distance_3d(coord2, coord4); l34 = distance_3d(coord3, coord4); Real s1, s2, s3, s4; s1 = (l12 + l23 + l13) * 0.5; s1 = sqrt(s1*(s1-l12)*(s1-l23)*(s1-l13)); s2 = (l12 + l24 + l14) * 0.5; s2 = sqrt(s2*(s2-l12)*(s2-l24)*(s2-l14)); s3 = (l23 + l34 + l24) * 0.5; s3 = sqrt(s3*(s3-l23)*(s3-l34)*(s3-l24)); s4 = (l13 + l34 + l24) * 0.5; s4 = sqrt(s4*(s4-l13)*(s4-l34)*(s4-l14)); Real volume = Math::tetrahedron_volume(coord1,coord2,coord3,coord4); return 3*volume/(s1+s2+s3+s4); } /* -------------------------------------------------------------------------- */ inline void Math::barycenter(const Real * coord, UInt nb_points, UInt spatial_dimension, Real * barycenter) { memset(barycenter, 0, spatial_dimension * sizeof(Real)); for (UInt n = 0; n < nb_points; ++n) { UInt offset = n * spatial_dimension; for (UInt i = 0; i < spatial_dimension; ++i) { barycenter[i] += coord[offset + i] / (Real) nb_points; } } } /* -------------------------------------------------------------------------- */ inline void Math::vector_2d(const Real * x, const Real * y, Real * res) { res[0] = y[0]-x[0]; res[1] = y[1]-x[1]; } /* -------------------------------------------------------------------------- */ inline void Math::vector_3d(const Real * x, const Real * y, Real * res) { res[0] = y[0]-x[0]; res[1] = y[1]-x[1]; res[2] = y[2]-x[2]; } +/* -------------------------------------------------------------------------- */ +inline bool Math::are_float_equal(const Real x, const Real y){ + return (fabs( x - y) < tolerance); +} diff --git a/fem/element_class.cc b/fem/element_class.cc index 666f8e111..eb84e66b2 100644 --- a/fem/element_class.cc +++ b/fem/element_class.cc @@ -1,157 +1,162 @@ /** * @file element_class.cc * @author Nicolas Richart * @date Tue Jul 20 10:12:44 2010 * * @brief Common part of element_classes * * @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.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template UInt ElementClass::nb_nodes_per_element = 0; template ElementType ElementClass::p1_element_type = _not_defined; template UInt ElementClass::nb_quadrature_points = 0; template UInt ElementClass::spatial_dimension = 0; template ElementType ElementClass::facet_type = _not_defined; + /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_point>::nb_nodes_per_element = 1; template<> ElementType ElementClass<_point>::p1_element_type = _point; +template<> UInt ElementClass<_point>::nb_quadrature_points = 1; +template<> Real ElementClass<_point>::quad[] = {0}; template<> ElementType ElementClass<_point>::facet_type = _not_defined; template<> UInt ElementClass<_point>::spatial_dimension = 0; template<> UInt ElementClass<_point>::nb_facets = 0; template<> UInt * ElementClass<_point>::facet_connectivity[] = {}; /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_segment_2>::nb_nodes_per_element = 2; template<> ElementType ElementClass<_segment_2>::p1_element_type = _segment_2; template<> UInt ElementClass<_segment_2>::nb_quadrature_points = 1; template<> Real ElementClass<_segment_2>::quad[] = {0}; template<> UInt ElementClass<_segment_2>::spatial_dimension = 1; template<> UInt ElementClass<_segment_2>::nb_facets = 2; template<> ElementType ElementClass<_segment_2>::facet_type = _point; template<> UInt ElementClass<_segment_2>::vec_facet_connectivity[]= {0, 1}; template<> UInt * ElementClass<_segment_2>::facet_connectivity[] = {&vec_facet_connectivity[0], &vec_facet_connectivity[1]}; /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_segment_3>::nb_nodes_per_element = 3; template<> ElementType ElementClass<_segment_3>::p1_element_type = _segment_2; template<> UInt ElementClass<_segment_3>::nb_quadrature_points = 2; template<> Real ElementClass<_segment_3>::quad[] = {-1./sqrt(3.), 1./sqrt(3.)}; template<> UInt ElementClass<_segment_3>::spatial_dimension = 1; template<> UInt ElementClass<_segment_3>::nb_facets = 2; template<> ElementType ElementClass<_segment_3>::facet_type = _point; template<> UInt ElementClass<_segment_3>::vec_facet_connectivity[]= {0, 1}; template<> UInt * ElementClass<_segment_3>::facet_connectivity[] = {&vec_facet_connectivity[0], &vec_facet_connectivity[1]}; /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_triangle_3>::nb_nodes_per_element = 3; template<> ElementType ElementClass<_triangle_3>::p1_element_type = _triangle_3; template<> UInt ElementClass<_triangle_3>::nb_quadrature_points = 1; template<> Real ElementClass<_triangle_3>::quad[] = {1./3., 1./3.}; template<> UInt ElementClass<_triangle_3>::spatial_dimension = 2; template<> UInt ElementClass<_triangle_3>::nb_facets = 3; template<> ElementType ElementClass<_triangle_3>::facet_type = _segment_2; template<> UInt ElementClass<_triangle_3>::vec_facet_connectivity[]= {0, 1, 1, 2, 2, 0}; template<> UInt * ElementClass<_triangle_3>::facet_connectivity[] = {&vec_facet_connectivity[0], &vec_facet_connectivity[2], &vec_facet_connectivity[4]}; /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_triangle_6>::nb_nodes_per_element = 6; template<> ElementType ElementClass<_triangle_6>::p1_element_type = _triangle_3; template<> UInt ElementClass<_triangle_6>::nb_quadrature_points = 3; template<> Real ElementClass<_triangle_6>::quad[] = {1./6., 1./6., 2./3., 1./6., 1./6., 2./3.}; template<> UInt ElementClass<_triangle_6>::spatial_dimension = 2; template<> UInt ElementClass<_triangle_6>::nb_facets = 3; template<> ElementType ElementClass<_triangle_6>::facet_type = _segment_3; template<> UInt ElementClass<_triangle_6>::vec_facet_connectivity[]= {0, 1, 3, 1, 2, 4, 2, 0, 5}; template<> UInt * ElementClass<_triangle_6>::facet_connectivity[] = {&vec_facet_connectivity[0], &vec_facet_connectivity[3], &vec_facet_connectivity[6]}; /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_tetrahedron_4>::nb_nodes_per_element = 4; template<> ElementType ElementClass<_tetrahedron_4>::p1_element_type = _tetrahedron_4; template<> UInt ElementClass<_tetrahedron_4>::nb_quadrature_points = 1; template<> Real ElementClass<_tetrahedron_4>::quad[] = {1./4., 1./4., 1./4.}; template<> UInt ElementClass<_tetrahedron_4>::spatial_dimension = 3; template<> UInt ElementClass<_tetrahedron_4>::nb_facets = 4; template<> ElementType ElementClass<_tetrahedron_4>::facet_type = _triangle_3; template<> UInt ElementClass<_tetrahedron_4>::vec_facet_connectivity[]= {0, 2, 1, 1, 2, 3, 2, 0, 3, 0, 1, 3}; template<> UInt * ElementClass<_tetrahedron_4>::facet_connectivity[] = {&vec_facet_connectivity[0], &vec_facet_connectivity[3], &vec_facet_connectivity[6], &vec_facet_connectivity[9]}; /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_tetrahedron_10>::nb_nodes_per_element = 10; template<> ElementType ElementClass<_tetrahedron_10>::p1_element_type = _tetrahedron_4; template<> UInt ElementClass<_tetrahedron_10>::nb_quadrature_points = 4; template<> Real ElementClass<_tetrahedron_10>::quad[] = {0.1381966011250, 0.1381966011250, 0.1381966011250, // a = (5-sqrt(5))/20 0.5854101966250, 0.1381966011250, 0.1381966011250, // b = (5+3*sqrt(5))/20 0.1381966011250, 0.5854101966250, 0.1381966011250, 0.1381966011250, 0.1381966011250, 0.5854101966250}; template<> UInt ElementClass<_tetrahedron_10>::spatial_dimension = 3; template<> UInt ElementClass<_tetrahedron_10>::nb_facets = 4; template<> ElementType ElementClass<_tetrahedron_10>::facet_type = _triangle_6; template<> UInt ElementClass<_tetrahedron_10>::vec_facet_connectivity[]= {0, 2, 1, 6, 5, 4, 1, 2, 3, 5, 9, 8, 2, 0, 3, 6, 7, 9, 0, 1, 3, 4, 8, 7}; template<> UInt * ElementClass<_tetrahedron_10>::facet_connectivity[] = {&vec_facet_connectivity[0], &vec_facet_connectivity[6], &vec_facet_connectivity[12], &vec_facet_connectivity[18]}; /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_quadrangle_4>::nb_nodes_per_element = 4; template<> ElementType ElementClass<_quadrangle_4>::p1_element_type = _quadrangle_4; template<> UInt ElementClass<_quadrangle_4>::nb_quadrature_points = 1; template<> Real ElementClass<_quadrangle_4>::quad[] = {0, 0}; template<> UInt ElementClass<_quadrangle_4>::spatial_dimension = 2; template<> UInt ElementClass<_quadrangle_4>::nb_facets = 4; template<> ElementType ElementClass<_quadrangle_4>::facet_type = _segment_2; template<> UInt ElementClass<_quadrangle_4>::vec_facet_connectivity[]= {0, 1, 1, 2, 2, 3, 3, 0}; template<> UInt * ElementClass<_quadrangle_4>::facet_connectivity[] = {&vec_facet_connectivity[0], &vec_facet_connectivity[2], &vec_facet_connectivity[4], &vec_facet_connectivity[6]}; /* -------------------------------------------------------------------------- */ + + __END_AKANTU__ diff --git a/fem/fem.cc b/fem/fem.cc index 240df4611..d85667e73 100644 --- a/fem/fem.cc +++ b/fem/fem.cc @@ -1,677 +1,655 @@ /** * @file fem.cc * @author Nicolas Richart + * @author Guillaume Anciaux * @date Fri Jul 16 11:03:02 2010 * * @brief Implementation of the FEM 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 "fem.hh" #include "mesh.hh" #include "element_class.hh" #include "aka_math.hh" + /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ + /* -------------------------------------------------------------------------- */ FEM::FEM(Mesh & mesh, UInt element_dimension, FEMID id, MemoryID memory_id) : Memory(memory_id), id(id) { AKANTU_DEBUG_IN(); this->element_dimension = (element_dimension != 0) ? element_dimension : mesh.getSpatialDimension(); init(); this->mesh = &mesh; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FEM::init() { for(UInt t = _not_defined; t < _max_element_type; ++t) { this->shapes [t] = NULL; this->shapes_derivatives[t] = NULL; this->jacobians [t] = NULL; this->normals_on_quad_points[t] = NULL; this->ghost_shapes [t] = NULL; this->ghost_shapes_derivatives[t] = NULL; this->ghost_jacobians [t] = NULL; } } /* -------------------------------------------------------------------------- */ FEM::~FEM() { AKANTU_DEBUG_IN(); mesh = NULL; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FEM::initShapeFunctions(GhostType ghost_type) { AKANTU_DEBUG_IN(); + Real * coord = mesh->getNodes().values; UInt spatial_dimension = mesh->getSpatialDimension(); const Mesh::ConnectivityTypeList & type_list = mesh->getConnectivityTypeList(ghost_type); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; UInt element_type_spatial_dimension = Mesh::getSpatialDimension(type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt size_of_shapes = FEM::getShapeSize(type); UInt size_of_shapesd = FEM::getShapeDerivativesSize(type); UInt nb_quadrature_points = FEM::getNbQuadraturePoints(type); if(element_type_spatial_dimension != element_dimension) continue; UInt * elem_val; UInt nb_element; std::string ghost = ""; if(ghost_type == _not_ghost) { elem_val = mesh->getConnectivity(type).values; nb_element = mesh->getConnectivity(type).getSize(); } else { ghost = "ghost_"; elem_val = mesh->getGhostConnectivity(type).values; nb_element = mesh->getGhostConnectivity(type).getSize(); } std::stringstream sstr_shapes; sstr_shapes << id << ":" << ghost << "shapes:" << type; Vector * shapes_tmp = &(alloc(sstr_shapes.str(), nb_element*nb_quadrature_points, size_of_shapes)); std::stringstream sstr_shapesd; sstr_shapesd << id << ":" << ghost << "shapes_derivatives:" << type; Vector * shapes_derivatives_tmp = &(alloc(sstr_shapesd.str(), nb_element*nb_quadrature_points, size_of_shapesd)); std::stringstream sstr_jacobians; sstr_jacobians << id << ":" << ghost << "jacobians:" << type; Vector * jacobians_tmp = &(alloc(sstr_jacobians.str(), nb_element*nb_quadrature_points, 1)); Real * shapes_val = shapes_tmp->values; Real * shapesd_val = shapes_derivatives_tmp->values; Real * jacobians_val = jacobians_tmp->values; - /* -------------------------------------------------------------------------- */ - /* compute shapes when no rotation is required */ - + /* ---------------------------------------------------------------------- */ #define COMPUTE_SHAPES(type) \ do { \ Real local_coord[spatial_dimension * nb_nodes_per_element]; \ for (UInt elem = 0; elem < nb_element; ++elem) { \ int offset = elem * nb_nodes_per_element; \ for (UInt id = 0; id < nb_nodes_per_element; ++id) { \ memcpy(local_coord + id * spatial_dimension, \ coord + elem_val[offset + id] * spatial_dimension, \ spatial_dimension*sizeof(Real)); \ } \ ElementClass::preComputeStandards(local_coord, \ spatial_dimension, \ shapes_val, \ shapesd_val, \ jacobians_val); \ shapes_val += size_of_shapes*nb_quadrature_points; \ shapesd_val += size_of_shapesd*nb_quadrature_points; \ jacobians_val += nb_quadrature_points; \ } \ } while(0) - -/* -------------------------------------------------------------------------- */ - - switch(type) { - case _segment_2 : { COMPUTE_SHAPES(_segment_2 ); break; } - case _segment_3 : { COMPUTE_SHAPES(_segment_3 ); break; } - case _triangle_3 : { COMPUTE_SHAPES(_triangle_3 ); break; } - case _triangle_6 : { COMPUTE_SHAPES(_triangle_6 ); break; } - case _tetrahedron_4 : { COMPUTE_SHAPES(_tetrahedron_4 ); break; } - case _tetrahedron_10 : { COMPUTE_SHAPES(_tetrahedron_10); break; } - case _quadrangle_4 : { COMPUTE_SHAPES(_quadrangle_4 ); break; } - case _point: - case _not_defined: - case _max_element_type: { - AKANTU_DEBUG_ERROR("Wrong type : " << type); - break; } - } + /* ---------------------------------------------------------------------- */ + + AKANTU_BOOST_ELEMENT_SWITCH(COMPUTE_SHAPES) #undef COMPUTE_SHAPES if(ghost_type == _not_ghost) { shapes[type] = shapes_tmp; shapes_derivatives[type] = shapes_derivatives_tmp; jacobians[type] = jacobians_tmp; } else { ghost_shapes[type] = shapes_tmp; ghost_shapes_derivatives[type] = shapes_derivatives_tmp; ghost_jacobians[type] = jacobians_tmp; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FEM::computeNormalsOnQuadPoints(GhostType ghost_type) { AKANTU_DEBUG_IN(); + Real * coord = mesh->getNodes().values; UInt spatial_dimension = mesh->getSpatialDimension(); + //allocate the normal arrays + if (ghost_type == _not_ghost) + mesh->initByElementTypeRealVector(normals_on_quad_points,spatial_dimension,element_dimension, + id,"normals_onquad",ghost_type); + else{ + AKANTU_DEBUG_ERROR("to be implemented"); + } + + //loop over the type to build the normals const Mesh::ConnectivityTypeList & type_list = mesh->getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; UInt element_type_spatial_dimension = Mesh::getSpatialDimension(type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quad_points = FEM::getNbQuadraturePoints(type); if(element_type_spatial_dimension != element_dimension) continue; UInt * elem_val; UInt nb_element; std::string ghost = ""; + Real * normals_on_quad_val = NULL; + if(ghost_type == _not_ghost) { elem_val = mesh->getConnectivity(type).values; nb_element = mesh->getConnectivity(type).getSize(); + normals_on_quad_val = normals_on_quad_points[type]->values; } else { ghost = "ghost_"; elem_val = mesh->getGhostConnectivity(type).values; nb_element = mesh->getGhostConnectivity(type).getSize(); } - std::stringstream sstr_normals_on_quad; - sstr_normals_on_quad << id << ":" << ghost << "normals_onquad:" << type; - Vector * normals_on_quad_tmp = &(alloc(sstr_normals_on_quad.str(), - nb_element*nb_quad_points, - spatial_dimension)); - - Real * normals_on_quad_val = normals_on_quad_tmp->values; + /* ---------------------------------------------------------------------- */ #define COMPUTE_NORMALS_ON_QUAD(type) \ do { \ Real local_coord[spatial_dimension * nb_nodes_per_element]; \ for (UInt elem = 0; elem < nb_element; ++elem) { \ int offset = elem * nb_nodes_per_element; \ for (UInt id = 0; id < nb_nodes_per_element; ++id) { \ memcpy(local_coord + id * spatial_dimension, \ coord + elem_val[offset + id] * spatial_dimension, \ spatial_dimension*sizeof(Real)); \ } \ ElementClass::computeNormalsOnQuadPoint(local_coord, \ spatial_dimension, \ normals_on_quad_val); \ normals_on_quad_val += spatial_dimension*nb_quad_points; \ } \ } while(0) - - switch(type) { - case _segment_2 : { COMPUTE_NORMALS_ON_QUAD(_segment_2 ); break; } - case _segment_3 : { COMPUTE_NORMALS_ON_QUAD(_segment_3 ); break; } - case _triangle_3 : { COMPUTE_NORMALS_ON_QUAD(_triangle_3 ); break; } - case _triangle_6 : { COMPUTE_NORMALS_ON_QUAD(_triangle_6 ); break; } - case _tetrahedron_4 : { COMPUTE_NORMALS_ON_QUAD(_tetrahedron_4 ); break; } - case _tetrahedron_10 : { COMPUTE_NORMALS_ON_QUAD(_tetrahedron_10); break; } - case _quadrangle_4 : { COMPUTE_NORMALS_ON_QUAD(_quadrangle_4 ); break; } - case _point: - case _not_defined: - case _max_element_type: { - AKANTU_DEBUG_ERROR("Wrong type : " << type); - break; } - } + /* ---------------------------------------------------------------------- */ + + AKANTU_BOOST_ELEMENT_SWITCH(COMPUTE_NORMALS_ON_QUAD) #undef COMPUTE_NORMALS_ON_QUAD - if(ghost_type == _not_ghost) { - normals_on_quad_points[type] = normals_on_quad_tmp; - } else { - AKANTU_DEBUG_ERROR("to be implemented"); - } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FEM::interpolateOnQuadraturePoints(const Vector &in_u, Vector &out_uq, UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type, const Vector * filter_elements) const { AKANTU_DEBUG_IN(); Vector * shapes_loc; UInt nb_element; UInt * conn_val; if(ghost_type == _not_ghost) { shapes_loc = shapes[type]; nb_element = mesh->getNbElement(type); conn_val = mesh->getConnectivity(type).values; } else { shapes_loc = ghost_shapes[type]; nb_element = mesh->getNbGhostElement(type); conn_val = mesh->getGhostConnectivity(type).values; } AKANTU_DEBUG_ASSERT(shapes_loc != NULL, "No shapes for the type " << type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = FEM::getNbQuadraturePoints(type); UInt size_of_shapes = FEM::getShapeSize(type); AKANTU_DEBUG_ASSERT(in_u.getSize() == mesh->getNbNodes(), "The vector in_u(" << in_u.getID() << ") has not the good size."); AKANTU_DEBUG_ASSERT(in_u.getNbComponent() == nb_degre_of_freedom, "The vector in_u(" << in_u.getID() << ") has not the good number of component."); AKANTU_DEBUG_ASSERT(out_uq.getNbComponent() == nb_degre_of_freedom , "The vector out_uq(" << out_uq.getID() << ") has not the good number of component."); UInt * filter_elem_val = NULL; if(filter_elements != NULL) { nb_element = filter_elements->getSize(); filter_elem_val = filter_elements->values; } out_uq.resize(nb_element * nb_quadrature_points); Real * shape_val = shapes_loc->values; Real * u_val = in_u.values; Real * uq_val = out_uq.values; UInt offset_uq = out_uq.getNbComponent()*nb_quadrature_points; Real * shape = shape_val; Real * u = static_cast(calloc(nb_nodes_per_element * nb_degre_of_freedom, sizeof(Real))); for (UInt el = 0; el < nb_element; ++el) { UInt el_offset = el * nb_nodes_per_element; if(filter_elements != NULL) { shape = shape_val + filter_elem_val[el] * size_of_shapes*nb_quadrature_points; el_offset = filter_elem_val[el] * nb_nodes_per_element; } for (UInt n = 0; n < nb_nodes_per_element; ++n) { memcpy(u + n * nb_degre_of_freedom, u_val + conn_val[el_offset + n] * nb_degre_of_freedom, nb_degre_of_freedom * sizeof(Real)); } /// Uq = Shape * U : matrix product Math::matrix_matrix(nb_quadrature_points, nb_degre_of_freedom, nb_nodes_per_element, shape, u, uq_val); uq_val += offset_uq; if(filter_elements == NULL) { shape += size_of_shapes*nb_quadrature_points; } } free(u); #undef INIT_VARIABLES AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FEM::gradientOnQuadraturePoints(const Vector &in_u, Vector &out_nablauq, UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type, const Vector * filter_elements) const { AKANTU_DEBUG_IN(); Vector * shapesd_loc; UInt nb_element; UInt * conn_val; if(ghost_type == _not_ghost) { shapesd_loc = shapes_derivatives[type]; nb_element = mesh->getNbElement(type); conn_val = mesh->getConnectivity(type).values; } else { shapesd_loc = ghost_shapes_derivatives[type]; nb_element = mesh->getNbGhostElement(type); conn_val = mesh->getGhostConnectivity(type).values; } AKANTU_DEBUG_ASSERT(shapesd_loc != NULL, "No shapes for the type " << type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt size_of_shapes_derivatives = FEM::getShapeDerivativesSize(type); UInt nb_quadrature_points = FEM::getNbQuadraturePoints(type); UInt * filter_elem_val = NULL; if(filter_elements != NULL) { nb_element = filter_elements->getSize(); filter_elem_val = filter_elements->values; } AKANTU_DEBUG_ASSERT(in_u.getSize() == mesh->getNbNodes(), "The vector in_u(" << in_u.getID() << ") has not the good size."); AKANTU_DEBUG_ASSERT(in_u.getNbComponent() == nb_degre_of_freedom , "The vector in_u(" << in_u.getID() << ") has not the good number of component."); AKANTU_DEBUG_ASSERT(out_nablauq.getNbComponent() == nb_degre_of_freedom * element_dimension, "The vector out_nablauq(" << out_nablauq.getID() << ") has not the good number of component."); out_nablauq.resize(nb_element * nb_quadrature_points); Real * shaped_val = shapesd_loc->values; Real * u_val = in_u.values; Real * nablauq_val = out_nablauq.values; UInt offset_nablauq = nb_degre_of_freedom * element_dimension; UInt offset_shaped = nb_nodes_per_element * element_dimension; Real * shaped = shaped_val; Real * u = static_cast(calloc(nb_nodes_per_element * nb_degre_of_freedom, sizeof(Real))); for (UInt el = 0; el < nb_element; ++el) { UInt el_offset = el * nb_nodes_per_element; if(filter_elements != NULL) { shaped = shaped_val + filter_elem_val[el] * size_of_shapes_derivatives*nb_quadrature_points; el_offset = filter_elem_val[el] * nb_nodes_per_element; } for (UInt n = 0; n < nb_nodes_per_element; ++n) { memcpy(u + n * nb_degre_of_freedom, u_val + conn_val[el_offset + n] * nb_degre_of_freedom, nb_degre_of_freedom * sizeof(Real)); } for (UInt q = 0; q < nb_quadrature_points; ++q) { /// \nabla(U) = U^t * dphi/dx Math::matrixt_matrix(nb_degre_of_freedom, element_dimension, nb_nodes_per_element, u, shaped, nablauq_val); nablauq_val += offset_nablauq; shaped += offset_shaped; } } free(u); #undef INIT_VARIABLES AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FEM::integrate(const Vector & in_f, Vector &intf, UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type, const Vector * filter_elements) const { AKANTU_DEBUG_IN(); Vector * jac_loc; UInt nb_element; if(ghost_type == _not_ghost) { jac_loc = jacobians[type]; nb_element = mesh->getNbElement(type); } else { jac_loc = ghost_jacobians[type]; nb_element = mesh->getNbGhostElement(type); } UInt nb_quadrature_points = FEM::getNbQuadraturePoints(type); UInt * filter_elem_val = NULL; if(filter_elements != NULL) { nb_element = filter_elements->getSize(); filter_elem_val = filter_elements->values; } AKANTU_DEBUG_ASSERT(in_f.getSize() == nb_element * nb_quadrature_points, "The vector in_f(" << in_f.getID() << " size " << in_f.getSize() << ") has not the good size (" << nb_element << ")."); AKANTU_DEBUG_ASSERT(in_f.getNbComponent() == nb_degre_of_freedom , "The vector in_f(" << in_f.getID() << ") has not the good number of component."); AKANTU_DEBUG_ASSERT(intf.getNbComponent() == nb_degre_of_freedom, "The vector intf(" << intf.getID() << ") has not the good number of component."); intf.resize(nb_element); Real * in_f_val = in_f.values; Real * intf_val = intf.values; Real * jac_val = jac_loc->values; UInt offset_in_f = in_f.getNbComponent()*nb_quadrature_points; UInt offset_intf = intf.getNbComponent(); Real * jac = jac_val; for (UInt el = 0; el < nb_element; ++el) { if(filter_elements != NULL) { jac = jac_val + filter_elem_val[el] * nb_quadrature_points; } integrate(in_f_val, jac, intf_val, nb_degre_of_freedom, nb_quadrature_points); in_f_val += offset_in_f; intf_val += offset_intf; if(filter_elements == NULL) { jac += nb_quadrature_points; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real FEM::integrate(const Vector & in_f, const ElementType & type, GhostType ghost_type, const Vector * filter_elements) const { AKANTU_DEBUG_IN(); Vector * jac_loc; UInt nb_element; if(ghost_type == _not_ghost) { jac_loc = jacobians[type]; nb_element = mesh->getNbElement(type); } else { jac_loc = ghost_jacobians[type]; nb_element = mesh->getNbGhostElement(type); } UInt nb_quadrature_points = FEM::getNbQuadraturePoints(type); UInt * filter_elem_val = NULL; if(filter_elements != NULL) { nb_element = filter_elements->getSize(); filter_elem_val = filter_elements->values; } AKANTU_DEBUG_ASSERT(in_f.getSize() == nb_element * nb_quadrature_points, "The vector in_f(" << in_f.getID() << ") has not the good size."); AKANTU_DEBUG_ASSERT(in_f.getNbComponent() == 1, "The vector in_f(" << in_f.getID() << ") has not the good number of component."); Real intf = 0.; Real * in_f_val = in_f.values; Real * jac_val = jac_loc->values; UInt offset_in_f = in_f.getNbComponent() * nb_quadrature_points; Real * jac = jac_val; for (UInt el = 0; el < nb_element; ++el) { if(filter_elements != NULL) { jac = jac_val + filter_elem_val[el] * nb_quadrature_points; } Real el_intf = 0; integrate(in_f_val, jac, &el_intf, 1, nb_quadrature_points); intf += el_intf; in_f_val += offset_in_f; if(filter_elements == NULL) { jac += nb_quadrature_points; } } AKANTU_DEBUG_OUT(); return intf; } /* -------------------------------------------------------------------------- */ void FEM::assembleVector(const Vector & elementary_vect, Vector & nodal_values, UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type, const Vector * filter_elements, Real scale_factor) const { AKANTU_DEBUG_IN(); UInt nb_element; UInt * conn_val; if(ghost_type == _not_ghost) { nb_element = mesh->getNbElement(type); conn_val = mesh->getConnectivity(type).values; } else { nb_element = mesh->getNbGhostElement(type); conn_val = mesh->getGhostConnectivity(type).values; } UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_nodes = mesh->getNbNodes(); UInt * filter_elem_val = NULL; if(filter_elements != NULL) { nb_element = filter_elements->getSize(); filter_elem_val = filter_elements->values; } AKANTU_DEBUG_ASSERT(elementary_vect.getSize() == nb_element, "The vector elementary_vect(" << elementary_vect.getID() << ") has not the good size."); AKANTU_DEBUG_ASSERT(elementary_vect.getNbComponent() == nb_degre_of_freedom*nb_nodes_per_element, "The vector elementary_vect(" << elementary_vect.getID() << ") has not the good number of component."); AKANTU_DEBUG_ASSERT(nodal_values.getNbComponent() == nb_degre_of_freedom, "The vector nodal_values(" << nodal_values.getID() << ") has not the good number of component."); nodal_values.resize(nb_nodes); Real * elementary_vect_val = elementary_vect.values; Real * nodal_values_val = nodal_values.values; for (UInt el = 0; el < nb_element; ++el) { UInt el_offset = el * nb_nodes_per_element; if(filter_elements != NULL) { el_offset = filter_elem_val[el] * nb_nodes_per_element; } for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = conn_val[el_offset + n]; UInt offset_node = node * nb_degre_of_freedom; for (UInt d = 0; d < nb_degre_of_freedom; ++d) { nodal_values_val[offset_node + d] += scale_factor * elementary_vect_val[d]; } elementary_vect_val += nb_degre_of_freedom; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FEM::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "FEM [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + element dimension : " << element_dimension << std::endl; stream << space << " + mesh [" << std::endl; mesh->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + connectivity type information [" << std::endl; const Mesh::ConnectivityTypeList & type_list = mesh->getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if (mesh->getSpatialDimension(*it) != element_dimension) continue; stream << space << AKANTU_INDENT << AKANTU_INDENT << " + " << *it <<" [" << std::endl; if(shapes[*it]) { shapes [*it]->printself(stream, indent + 3); shapes_derivatives[*it]->printself(stream, indent + 3); jacobians [*it]->printself(stream, indent + 3); } stream << space << AKANTU_INDENT << AKANTU_INDENT << "]" << std::endl; } stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/fem/fem.hh b/fem/fem.hh index 41047e1b9..7cf87b981 100644 --- a/fem/fem.hh +++ b/fem/fem.hh @@ -1,222 +1,223 @@ /** * @file fem.hh * @author Nicolas Richart * @date Fri Jul 16 10:24:24 2010 * * @brief FEM 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_FEM_HH__ #define __AKANTU_FEM_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "mesh.hh" #include "element_class.hh" -/* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ + + + class FEM : public Memory { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: FEM(Mesh & mesh, UInt spatial_dimension = 0, FEMID id = "fem", MemoryID memory_id = 0); virtual ~FEM(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// pre-compute all the shape functions, their derivatives and the jacobians void initShapeFunctions(GhostType ghost_type = _not_ghost); /// build the profile of the sparse matrix corresponding to the mesh void initSparseMatrixProfile(SparseMatrixType sparse_matrix_type = _unsymmetric); /// pre-compute normals on quadrature points void computeNormalsOnQuadPoints(GhostType ghost_type = _not_ghost); /// interpolate nodal values on the quadrature points void interpolateOnQuadraturePoints(const Vector &u, Vector &uq, UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL) const; /// compute the gradient of u on the quadrature points void gradientOnQuadraturePoints(const Vector &u, Vector &nablauq, UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL) const; /// integrate f for all elements of type "type" void integrate(const Vector & f, Vector &intf, UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL) const; /// integrate f on the element "elem" of type "type" inline void integrate(const Vector & f, Real *intf, UInt nb_degre_of_freedom, const Element & elem, GhostType ghost_type = _not_ghost) const; /// integrate a scalar value on all elements of type "type" Real integrate(const Vector & f, const ElementType & type, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL) const; /// assemble vectors void assembleVector(const Vector & elementary_vect, Vector & nodal_values, UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL, Real scale_factor = 1) const; /// assemble matrix in the complete sparse matrix void assembleMatrix() {}; /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; private: /// initialise the class void init(); /// integrate on one element inline void integrate(Real *f, Real *jac, Real * inte, UInt nb_degre_of_freedom, UInt nb_quadrature_points) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(ElementDimension, element_dimension, UInt); /// get the mesh contained in the fem object inline Mesh & getMesh() const; /// get the size of the shapes returned by the element class static inline UInt getShapeSize(const ElementType & type); /// get the number of quadrature points static inline UInt getNbQuadraturePoints(const ElementType & type); /// get the size of the shapes derivatives returned by the element class static inline UInt getShapeDerivativesSize(const ElementType & type); /// get the in-radius of an element static inline Real getElementInradius(Real * coord, const ElementType & type); /// get a the shapes vector AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Shapes, shapes, const Vector &); /// get a the shapes derivatives vector AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ShapesDerivatives, shapes_derivatives, const Vector &); /// get the normals on quadrature points AKANTU_GET_MACRO_BY_ELEMENT_TYPE(NormalsOnQuadPoints, normals_on_quad_points, const Vector &); /// get a the ghost shapes vector AKANTU_GET_MACRO_BY_ELEMENT_TYPE(GhostShapes, ghost_shapes, const Vector &); /// get a the ghost shapes derivatives vector AKANTU_GET_MACRO_BY_ELEMENT_TYPE(GhostShapesDerivatives, ghost_shapes_derivatives, const Vector &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// id of the fem object FEMID id; /// spatial dimension of the problem UInt element_dimension; /// the mesh on which all computation are made Mesh * mesh; /// shape functions for all elements ByElementTypeReal shapes; /// shape derivatives for all elements ByElementTypeReal shapes_derivatives; /// jacobians for all elements ByElementTypeReal jacobians; /// normals at quadrature points ByElementTypeReal normals_on_quad_points; /// shape functions for all elements ByElementTypeReal ghost_shapes; /// shape derivatives for all elements ByElementTypeReal ghost_shapes_derivatives; /// jacobians for all elements ByElementTypeReal ghost_jacobians; }; - /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "fem_inline_impl.cc" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const FEM & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_FEM_HH__ */ diff --git a/fem/fem_inline_impl.cc b/fem/fem_inline_impl.cc index 285096f3d..95252c9ff 100644 --- a/fem/fem_inline_impl.cc +++ b/fem/fem_inline_impl.cc @@ -1,192 +1,137 @@ /** * @file fem_inline_impl.cc * @author Nicolas Richart + * @author Guillaume Anciaux * @date Mon Jul 19 12:21:36 2010 * * @brief Implementation of the inline functions of the FEM 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 . * */ /* -------------------------------------------------------------------------- */ inline void FEM::integrate(const Vector & f, Real * intf, UInt nb_degre_of_freedom, const Element & elem, GhostType ghost_type) const { Vector * jac_loc; if(ghost_type == _not_ghost) { jac_loc = jacobians[elem.type]; } else { jac_loc = ghost_jacobians[elem.type]; } UInt nb_quadrature_points = FEM::getNbQuadraturePoints(elem.type); AKANTU_DEBUG_ASSERT(f.getNbComponent() == nb_degre_of_freedom , "The vector f do not have the good number of component."); Real * f_val = f.values + elem.element * f.getNbComponent(); Real * jac_val = jac_loc->values + elem.element * nb_quadrature_points; integrate(f_val, jac_val, intf, nb_degre_of_freedom, nb_quadrature_points); } /* -------------------------------------------------------------------------- */ inline void FEM::integrate(Real *f, Real *jac, Real * inte, UInt nb_degre_of_freedom, UInt nb_quadrature_points) const { memset(inte, 0, nb_degre_of_freedom * sizeof(Real)); Real *cjac = jac; for (UInt q = 0; q < nb_quadrature_points; ++q) { for (UInt dof = 0; dof < nb_degre_of_freedom; ++dof) { inte[dof] += *f * *cjac; ++f; } ++cjac; } } /* -------------------------------------------------------------------------- */ inline Mesh & FEM::getMesh() const { return *mesh; } /* -------------------------------------------------------------------------- */ inline UInt FEM::getNbQuadraturePoints(const ElementType & type) { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = 0; #define GET_NB_QUAD_POINTS(type) \ nb_quadrature_points = ElementClass::getNbQuadraturePoints() - switch(type) { - case _segment_2 : { GET_NB_QUAD_POINTS(_segment_2 ); break; } - case _segment_3 : { GET_NB_QUAD_POINTS(_segment_3 ); break; } - case _triangle_3 : { GET_NB_QUAD_POINTS(_triangle_3 ); break; } - case _triangle_6 : { GET_NB_QUAD_POINTS(_triangle_6 ); break; } - case _tetrahedron_4 : { GET_NB_QUAD_POINTS(_tetrahedron_4 ); break; } - case _tetrahedron_10 : { GET_NB_QUAD_POINTS(_tetrahedron_10); break; } - case _quadrangle_4 : { GET_NB_QUAD_POINTS(_quadrangle_4 ); break; } - case _point: - case _not_defined: - case _max_element_type: { - AKANTU_DEBUG_ERROR("Wrong type : " << type); - break; } - } - + AKANTU_BOOST_ELEMENT_SWITCH(GET_NB_QUAD_POINTS) #undef GET_NB_QUAD_POINTS AKANTU_DEBUG_OUT(); return nb_quadrature_points; } /* -------------------------------------------------------------------------- */ inline UInt FEM::getShapeSize(const ElementType & type) { AKANTU_DEBUG_IN(); UInt shape_size = 0; #define GET_SHAPE_SIZE(type) \ shape_size = ElementClass::getShapeSize() - switch(type) { - case _segment_2 : { GET_SHAPE_SIZE(_segment_2 ); break; } - case _segment_3 : { GET_SHAPE_SIZE(_segment_3 ); break; } - case _triangle_3 : { GET_SHAPE_SIZE(_triangle_3 ); break; } - case _triangle_6 : { GET_SHAPE_SIZE(_triangle_6 ); break; } - case _tetrahedron_4 : { GET_SHAPE_SIZE(_tetrahedron_4 ); break; } - case _tetrahedron_10 : { GET_SHAPE_SIZE(_tetrahedron_10); break; } - case _quadrangle_4 : { GET_SHAPE_SIZE(_quadrangle_4 ); break; } - case _point: - case _not_defined: - case _max_element_type: { - AKANTU_DEBUG_ERROR("Wrong type : " << type); - break; } - } - + AKANTU_BOOST_ELEMENT_SWITCH(GET_SHAPE_SIZE) #undef GET_SHAPE_SIZE AKANTU_DEBUG_OUT(); return shape_size; } /* -------------------------------------------------------------------------- */ inline UInt FEM::getShapeDerivativesSize(const ElementType & type) { AKANTU_DEBUG_IN(); UInt shape_derivatives_size = 0; #define GET_SHAPE_DERIVATIVES_SIZE(type) \ shape_derivatives_size = ElementClass::getShapeDerivativesSize() - switch(type) { - case _segment_2 : { GET_SHAPE_DERIVATIVES_SIZE(_segment_2 ); break; } - case _segment_3 : { GET_SHAPE_DERIVATIVES_SIZE(_segment_3 ); break; } - case _triangle_3 : { GET_SHAPE_DERIVATIVES_SIZE(_triangle_3 ); break; } - case _triangle_6 : { GET_SHAPE_DERIVATIVES_SIZE(_triangle_6 ); break; } - case _tetrahedron_4 : { GET_SHAPE_DERIVATIVES_SIZE(_tetrahedron_4 ); break; } - case _tetrahedron_10 : { GET_SHAPE_DERIVATIVES_SIZE(_tetrahedron_10); break; } - case _quadrangle_4 : { GET_SHAPE_DERIVATIVES_SIZE(_quadrangle_4 ); break; } - case _point: - case _not_defined: - case _max_element_type: { - AKANTU_DEBUG_ERROR("Wrong type : " << type); - break; } - } - + AKANTU_BOOST_ELEMENT_SWITCH(GET_SHAPE_DERIVATIVES_SIZE) #undef GET_SHAPE_DERIVATIVES_SIZE AKANTU_DEBUG_OUT(); return shape_derivatives_size; } /* -------------------------------------------------------------------------- */ inline Real FEM::getElementInradius(Real * coord, const ElementType & type) { AKANTU_DEBUG_IN(); Real inradius = 0; #define GET_INRADIUS(type) \ inradius = ElementClass::getInradius(coord); \ - switch(type) { - case _segment_2 : { GET_INRADIUS(_segment_2 ); break; } - case _segment_3 : { GET_INRADIUS(_segment_3 ); break; } - case _triangle_3 : { GET_INRADIUS(_triangle_3 ); break; } - case _triangle_6 : { GET_INRADIUS(_triangle_6 ); break; } - case _tetrahedron_4 : { GET_INRADIUS(_tetrahedron_4 ); break; } - case _tetrahedron_10 : { GET_INRADIUS(_tetrahedron_10); break; } - case _quadrangle_4 : { GET_INRADIUS(_quadrangle_4 ); break; } - case _point: - case _not_defined: - case _max_element_type: { - AKANTU_DEBUG_ERROR("Wrong type : " << type); - break; } - } - + AKANTU_BOOST_ELEMENT_SWITCH(GET_INRADIUS) #undef GET_INRADIUS AKANTU_DEBUG_OUT(); return inradius; } /* -------------------------------------------------------------------------- */ diff --git a/fem/mesh.cc b/fem/mesh.cc index 8f3a443a4..2c651e308 100644 --- a/fem/mesh.cc +++ b/fem/mesh.cc @@ -1,163 +1,216 @@ /** * @file mesh.cc * @author Nicolas Richart * @date Wed Jun 16 12:02:26 2010 * * @brief class handling meshes * * @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 /* -------------------------------------------------------------------------- */ #include "mesh.hh" #include "element_class.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ void Element::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "Element [" << type << ", " << element << "]"; } /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, const MeshID & id, const MemoryID & memory_id) : Memory(memory_id), id(id), nodes_global_ids(NULL), created_nodes(true), spatial_dimension(spatial_dimension), internal_facets_mesh(NULL), types_offsets(Vector(_max_element_type + 1, 1)), ghost_types_offsets(Vector(_max_element_type + 1, 1)), nb_surfaces(0) { AKANTU_DEBUG_IN(); initConnectivities(); std::stringstream sstr; sstr << id << ":coordinates"; this->nodes = &(alloc(sstr.str(), 0, this->spatial_dimension)); nb_global_nodes = 0; + computeBoundingBox(); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, const VectorID & nodes_id, const MeshID & id, const MemoryID & memory_id) : Memory(memory_id), id(id), nodes_global_ids(NULL), created_nodes(false), spatial_dimension(spatial_dimension), internal_facets_mesh(NULL), types_offsets(Vector(_max_element_type + 1, 1)), ghost_types_offsets(Vector(_max_element_type + 1, 1)) { AKANTU_DEBUG_IN(); initConnectivities(); this->nodes = &(getVector(nodes_id)); nb_global_nodes = nodes->getSize(); + computeBoundingBox(); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, Vector & nodes, const MeshID & id, const MemoryID & memory_id) : Memory(memory_id), id(id), created_nodes(false), spatial_dimension(spatial_dimension), internal_facets_mesh(NULL), types_offsets(Vector(_max_element_type + 1, 1)), ghost_types_offsets(Vector(_max_element_type + 1, 1)) { AKANTU_DEBUG_IN(); initConnectivities(); this->nodes = &(nodes); nb_global_nodes = nodes.getSize(); + computeBoundingBox(); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Mesh::initConnectivities() { for(UInt t = _not_defined; t < _max_element_type; ++t) { connectivities[t] = NULL; ghost_connectivities[t] = NULL; surface_id[t] = NULL; } this->types_offsets.resize(_max_element_type); } /* -------------------------------------------------------------------------- */ Mesh::~Mesh() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } -/* -------------------------------------------------------------------------- */ -// Vector & Mesh::createNormals(ElementType type) { -// AKANTU_DEBUG_IN(); -// AKANTU_DEBUG_ERROR("TOBEIMPLEMENTED"); -// AKANTU_DEBUG_OUT(); -// return *normals[type]; -// } - /* -------------------------------------------------------------------------- */ void Mesh::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "Mesh [" << std::endl; stream << space << " + id : " << this->id << std::endl; stream << space << " + spatial dimension : " << this->spatial_dimension << std::endl; stream << space << " + nodes [" << std::endl; nodes->printself(stream, indent+2); stream << space << " ]" << std::endl; ConnectivityTypeList::const_iterator it; for(it = type_set.begin(); it != type_set.end(); ++it) { stream << space << " + connectivities ("<< *it <<") [" << std::endl; (connectivities[*it])->printself(stream, indent+2); stream << space << " ]" << std::endl; } for(it = ghost_type_set.begin(); it != ghost_type_set.end(); ++it) { stream << space << " + ghost_connectivities ("<< *it <<") [" << std::endl; (ghost_connectivities[*it])->printself(stream, indent+2); stream << space << " ]" << std::endl; } stream << space << "]" << std::endl; } +/* -------------------------------------------------------------------------- */ + +void Mesh::computeBoundingBox(){ + AKANTU_DEBUG_IN(); + UInt dim = spatial_dimension; + memset(xmin,REAL_INIT_VALUE,3*sizeof(Real)); + memset(xmax,REAL_INIT_VALUE,3*sizeof(Real)); + + for (UInt k = 0; k < dim; ++k) { + xmin[k] = std::numeric_limits::max(); + xmax[k] = std::numeric_limits::min(); + } + + Real * coords = nodes->values; + for (UInt i = 0; i < nodes->getSize(); ++i) { + for (UInt k = 0; k < dim; ++k) { + xmin[k] = fmin(xmin[k],coords[dim*i+k]); + xmax[k] = fmax(xmax[k],coords[dim*i+k]); + } + } + AKANTU_DEBUG_OUT(); +} /* -------------------------------------------------------------------------- */ +void Mesh::initByElementTypeRealVector(ByElementTypeReal & vect, + UInt nb_component, + UInt dim, + const std::string & obj_id, + const std::string & vect_id, + GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + for(UInt t = _not_defined; t < _max_element_type; ++t) + vect[t] = NULL; + + std::string ghost_id = ""; + + if (ghost_type == _ghost) { + ghost_id = "ghost_"; + } + + const Mesh::ConnectivityTypeList & type_list = getConnectivityTypeList(); + Mesh::ConnectivityTypeList::const_iterator it; + for(it = type_list.begin(); it != type_list.end(); ++it) { + if(Mesh::getSpatialDimension(*it) != dim) continue; + std::stringstream sstr_damage; sstr_damage << obj_id << ":" << ghost_id << vect_id << ":" << *it; + if (vect[*it] == NULL){ + vect[*it] = &(alloc(sstr_damage.str(), 0, + nb_component, REAL_INIT_VALUE)); + } + } + + AKANTU_DEBUG_OUT(); +} + + + __END_AKANTU__ diff --git a/fem/mesh.hh b/fem/mesh.hh index ad100519c..9c61ba348 100644 --- a/fem/mesh.hh +++ b/fem/mesh.hh @@ -1,338 +1,347 @@ /** * @file mesh.hh * @author Nicolas Richart * @date Wed Jun 16 11:53:53 2010 * * @brief the class representing the meshes * * @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_MESH_HH__ #define __AKANTU_MESH_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "aka_vector.hh" #include "element_class.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /// to store data Vector by element type typedef Vector * ByElementTypeReal[_max_element_type]; /// to store data Vector by element type typedef Vector * ByElementTypeInt [_max_element_type]; /// to store data Vector by element type typedef Vector * ByElementTypeUInt[_max_element_type]; /* -------------------------------------------------------------------------- */ class Element { public: Element(ElementType type = _not_defined, UInt element = 0) : type(type), element(element) {}; Element(const Element & element) { this->type = element.type; this->element = element.element; } ~Element() {}; /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; public: ElementType type; UInt element; }; /* -------------------------------------------------------------------------- */ /** * @class Mesh this contain the coordinates of the nodes in the Mesh.nodes Vector, * and the connectivity. The connectivity are stored in by element types. * * To know all the element types present in a mesh you can get the Mesh::ConnectivityTypeList * * In order to loop on all element you have to loop on all types like this : * @code const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { UInt nb_element = mesh.getNbElement(*it); UInt * conn = mesh.getConnectivity(*it).values; for(UInt e = 0; e < nb_element; ++e) { ... } } @endcode */ class Mesh : protected Memory { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /// constructor that create nodes coordinates array Mesh(UInt spatial_dimension, const MeshID & id = "mesh", const MemoryID & memory_id = 0); /// constructor that use an existing nodes coordinates array, by knowing its ID Mesh(UInt spatial_dimension, const VectorID & nodes_id, const MeshID & id = "mesh", const MemoryID & memory_id = 0); /** * constructor that use an existing nodes coordinates * array, by getting the vector of coordinates */ Mesh(UInt spatial_dimension, Vector & nodes, const MeshID & id = "mesh", const MemoryID & memory_id = 0); virtual ~Mesh(); /// @typedef ConnectivityTypeList list of the types present in a Mesh typedef std::set ConnectivityTypeList; typedef Vector * NormalsMap[_max_element_type]; private: /// initialize the connectivity to NULL void initConnectivities(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; - // Vector & createConnectivity(ElementType type, UInt nb_element); - - // Vector & createGhostConnectivity(ElementType type, UInt nb_element); - - /// create the array of normals - // Vector & createNormals(ElementType type); - + /// function that computes the bounding box (fills xmin, xmax) + void computeBoundingBox(); + + /// init a by-element-type real vector with provided ids + void initByElementTypeRealVector(ByElementTypeReal & v,UInt nb_component, + UInt dimension, + const std::string & obj_id, + const std::string & vec_id, + GhostType ghost_type); + /// convert a element to a linearized element inline UInt elementToLinearized(const Element & elem); /// convert a linearized element to an element inline Element linearizedToElement (UInt linearized_element); /// update the types offsets array for the conversions inline void updateTypesOffsets(); /// convert a element to a linearized element inline UInt ghostElementToLinearized(const Element & elem); /// convert a linearized element to an element inline Element ghostLinearizedToElement (UInt linearized_element); /// update the types offsets array for the conversions inline void updateGhostTypesOffsets(); + + + /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(ID, id, const MeshID &); /// get the spatial dimension of the mesh = number of component of the coordinates AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the nodes Vector aka coordinates AKANTU_GET_MACRO(Nodes, *nodes, const Vector &); /// get the number of nodes AKANTU_GET_MACRO(NbNodes, nodes->getSize(), UInt); /// get the Vector of global ids of the nodes (only used in parallel) AKANTU_GET_MACRO(GlobalNodesIds, *nodes_global_ids, const Vector &); /// get the global number of nodes AKANTU_GET_MACRO(NbGlobalNodes, nb_global_nodes, UInt); /// get the number of surfaces AKANTU_GET_MACRO(NbSurfaces, nb_surfaces, UInt); /// get the connectivity Vector for a given type AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Connectivity, connectivities, const Vector &); /// get the connecticity of ghost elements of a given type AKANTU_GET_MACRO_BY_ELEMENT_TYPE(GhostConnectivity, ghost_connectivities, const Vector &); // AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Normals, normals, const Vector &); /// get the number of element of a type in the mesh inline UInt getNbElement(const ElementType & type) const; /// get the number of ghost element of a type in the mesh inline UInt getNbGhostElement(const ElementType & type) const; /// get the connectivity list either for the elements or the ghost elements inline const ConnectivityTypeList & getConnectivityTypeList(GhostType ghost_type = _not_ghost) const; /// get the mesh of the internal facets inline const Mesh & getInternalFacetsMesh() const; /// compute the barycenter of a given element inline void getBarycenter(UInt element, ElementType type, Real * barycenter, GhostType ghost_type = _not_ghost) const; /// get the surface values of facets inline const Vector & getSurfaceId(const ElementType & type) const; /* ------------------------------------------------------------------------ */ /* Wrappers on ElementClass functions */ /* ------------------------------------------------------------------------ */ public: /// get the number of nodes per element for a given element type static inline UInt getNbNodesPerElement(const ElementType & type); /// get the number of nodes per element for a given element type considered as /// a first order element static inline ElementType getP1ElementType(const ElementType & type); /// get spatial dimension of a type of element static inline UInt getSpatialDimension(const ElementType & type); /// get number of facets of a given element type static inline UInt getNbFacetsPerElement(const ElementType & type); /// get number of facets of a given element type static inline UInt ** getFacetLocalConnectivity(const ElementType & type); /// get the type of the surface element associated to a given element static inline ElementType getFacetElementType(const ElementType & type); private: friend class MeshIOMSH; friend class MeshUtils; friend class Communicator; AKANTU_GET_MACRO(NodesPointer, nodes, Vector *); inline Vector * getNodesGlobalIdsPointer(); inline Vector * getConnectivityPointer(ElementType type); inline Vector * getGhostConnectivityPointer(ElementType type); inline Mesh * getInternalFacetsMeshPointer(); // inline Vector * getNormalsPointer(ElementType type) const; inline Vector * getSurfaceIdPointer(ElementType type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// id of the mesh MeshID id; /// array of the nodes coordinates Vector * nodes; /// global node ids Vector * nodes_global_ids; /// global number of nodes; UInt nb_global_nodes; /// boolean to know if the nodes have to be deleted with the mesh or not bool created_nodes; /// all class of elements present in this mesh (for heterogenous meshes) ByElementTypeUInt connectivities; /// map to normals for all class of elements present in this mesh ByElementTypeReal normals; /// list of all existing types in the mesh ConnectivityTypeList type_set; /// the spatial dimension of this mesh UInt spatial_dimension; /// internal facets mesh Mesh * internal_facets_mesh; /// types offsets Vector types_offsets; /// all class of elements present in this mesh (for heterogenous meshes) ByElementTypeUInt ghost_connectivities; /// list of all existing types in the mesh ConnectivityTypeList ghost_type_set; /// ghost types offsets Vector ghost_types_offsets; /// number of surfaces present in this mesh UInt nb_surfaces; /// surface id of the surface elements in this mesh ByElementTypeUInt surface_id; -}; + /// min of coordinates + Real xmin[3]; + /// max of coordinates + Real xmax[3]; +}; /* -------------------------------------------------------------------------- */ /* Inline functions */ /* -------------------------------------------------------------------------- */ /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Element & _this) { _this.printself(stream); return stream; } - #include "mesh_inline_impl.cc" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Mesh & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_MESH_HH__ */ diff --git a/fem/mesh_inline_impl.cc b/fem/mesh_inline_impl.cc index dac061a88..c8b500d33 100644 --- a/fem/mesh_inline_impl.cc +++ b/fem/mesh_inline_impl.cc @@ -1,460 +1,368 @@ /** * @file mesh_inline_impl.cc * @author Nicolas Richart * @date Wed Jul 14 23:58:08 2010 * * @brief Implementation of the inline functions of the mesh 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 . * */ /* -------------------------------------------------------------------------- */ inline UInt Mesh::elementToLinearized(const Element & elem) { AKANTU_DEBUG_ASSERT(elem.type < _max_element_type && elem.element < types_offsets.values[elem.type+1], "The element " << elem << "does not exists in the mesh " << id); return types_offsets.values[elem.type] + elem.element; } /* -------------------------------------------------------------------------- */ inline Element Mesh::linearizedToElement (UInt linearized_element) { UInt t; for (t = _not_defined + 1; linearized_element > types_offsets.values[t] && t <= _max_element_type; ++t); AKANTU_DEBUG_ASSERT(t < _max_element_type, "The linearized element " << linearized_element << "does not exists in the mesh " << id); return Element((ElementType) t, linearized_element-types_offsets.values[t]); } /* -------------------------------------------------------------------------- */ inline void Mesh::updateTypesOffsets() { UInt count = 0; for (UInt t = _not_defined; t <= _max_element_type; ++t) { types_offsets.values[t] = count; count += (t == _max_element_type || connectivities[t] == NULL) ? 0 : connectivities[t]->getSize(); } } /* -------------------------------------------------------------------------- */ inline UInt Mesh::ghostElementToLinearized(const Element & elem) { AKANTU_DEBUG_ASSERT(elem.type < _max_element_type && elem.element < ghost_types_offsets.values[elem.type+1], "The ghost element " << elem << "does not exists in the mesh " << id); return ghost_types_offsets.values[elem.type] + elem.element + types_offsets.values[_max_element_type]; } /* -------------------------------------------------------------------------- */ inline Element Mesh::ghostLinearizedToElement (UInt linearized_element) { AKANTU_DEBUG_ASSERT(linearized_element >= types_offsets.values[_max_element_type], "The linearized element " << linearized_element << "is not a ghost element in the mesh " << id); linearized_element -= types_offsets.values[_max_element_type]; UInt t; for (t = _not_defined + 1; linearized_element > ghost_types_offsets.values[t] && t <= _max_element_type; ++t); AKANTU_DEBUG_ASSERT(t < _max_element_type, "The ghost linearized element " << linearized_element << "does not exists in the mesh " << id); t--; return Element((ElementType) t, linearized_element - ghost_types_offsets.values[t]); } /* -------------------------------------------------------------------------- */ inline void Mesh::updateGhostTypesOffsets() { UInt count = 0; for (UInt t = _not_defined; t <= _max_element_type; ++t) { ghost_types_offsets.values[t] = count; count += (t == _max_element_type || ghost_connectivities[t] == NULL) ? 0 : ghost_connectivities[t]->getSize(); } } /* -------------------------------------------------------------------------- */ inline const Mesh::ConnectivityTypeList & Mesh::getConnectivityTypeList(GhostType ghost_type) const { if(ghost_type == _not_ghost) return type_set; else return ghost_type_set; } /* -------------------------------------------------------------------------- */ inline Vector * Mesh::getNodesGlobalIdsPointer() { AKANTU_DEBUG_IN(); if(nodes_global_ids == NULL) { std::stringstream sstr; sstr << id << ":nodes_global_ids"; nodes_global_ids = &(alloc(sstr.str(), nodes->getSize(), 1)); } AKANTU_DEBUG_OUT(); return nodes_global_ids; } /* -------------------------------------------------------------------------- */ inline Vector * Mesh::getConnectivityPointer(ElementType type) { AKANTU_DEBUG_IN(); if(connectivities[type] == NULL) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); std::stringstream sstr; sstr << id << ":connectivity:" << type; connectivities[type] = &(alloc(sstr.str(), 0, nb_nodes_per_element)); type_set.insert(type); AKANTU_DEBUG_INFO("The connectivity vector for the type " << type << " created"); updateTypesOffsets(); } AKANTU_DEBUG_OUT(); return connectivities[type]; } /* -------------------------------------------------------------------------- */ inline Vector * Mesh::getGhostConnectivityPointer(ElementType type) { AKANTU_DEBUG_IN(); if(ghost_connectivities[type] == NULL) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); std::stringstream sstr; sstr << id << ":ghost_connectivity:" << type; ghost_connectivities[type] = &(alloc(sstr.str(), 0, nb_nodes_per_element)); ghost_type_set.insert(type); AKANTU_DEBUG_INFO("The connectivity vector for the type " << type << " created"); updateGhostTypesOffsets(); } AKANTU_DEBUG_OUT(); return ghost_connectivities[type]; } -/* -------------------------------------------------------------------------- */ -// inline Vector * Mesh::getNormalsPointer(ElementType type) const { -// AKANTU_DEBUG_IN(); - -// AKANTU_DEBUG_OUT(); -// return normals[type]; -// } - /* -------------------------------------------------------------------------- */ inline const Mesh & Mesh::getInternalFacetsMesh() const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); if (!internal_facets_mesh) AKANTU_DEBUG_ERROR("internal facets mesh was not created before access => use mesh utils to that purpose"); return *internal_facets_mesh; } /* -------------------------------------------------------------------------- */ inline Mesh * Mesh::getInternalFacetsMeshPointer() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); if (!internal_facets_mesh){ std::stringstream name(this->id); name << ":internalfacets"; internal_facets_mesh = new Mesh(this->spatial_dimension-1,*this->nodes,name.str()); } return internal_facets_mesh; } /* -------------------------------------------------------------------------- */ inline Vector * Mesh::getSurfaceIdPointer(ElementType type) { AKANTU_DEBUG_IN(); if(surface_id[type] == NULL) { std::stringstream sstr; sstr << id << ":surface_id:" << type; surface_id[type] = &(alloc(sstr.str(), 0, 1)); AKANTU_DEBUG_INFO("The surface id vector for the type " << type << " created"); } AKANTU_DEBUG_OUT(); return surface_id[type]; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbElement(const ElementType & type) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(connectivities[type] != NULL, "No element of kind : " << type << " in " << id); AKANTU_DEBUG_OUT(); return connectivities[type]->getSize(); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbGhostElement(const ElementType & type) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(ghost_connectivities[type] != NULL, "No element of kind : " << type << " in " << id); AKANTU_DEBUG_OUT(); return ghost_connectivities[type]->getSize(); } /* -------------------------------------------------------------------------- */ inline void Mesh::getBarycenter(UInt element, ElementType type, Real * barycenter, GhostType ghost_type) const { AKANTU_DEBUG_IN(); UInt * conn_val; if (ghost_type == _not_ghost) { conn_val = getConnectivity(type).values; } else { conn_val = getGhostConnectivity(type).values; } UInt nb_nodes_per_element = getNbNodesPerElement(type); Real local_coord[spatial_dimension * nb_nodes_per_element]; UInt offset = element * nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) { memcpy(local_coord + n * spatial_dimension, nodes->values + conn_val[offset + n] * spatial_dimension, spatial_dimension*sizeof(Real)); } Math::barycenter(local_coord, nb_nodes_per_element, spatial_dimension, barycenter); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline const Vector & Mesh::getSurfaceId(const ElementType & type) const{ AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(surface_id[type] != NULL, "No element of kind : " << type << " in " << id); AKANTU_DEBUG_OUT(); return *surface_id[type]; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbNodesPerElement(const ElementType & type) { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = 0; #define GET_NB_NODES_PER_ELEMENT(type) \ nb_nodes_per_element = ElementClass::getNbNodesPerElement() - switch(type) { - case _segment_2 : { GET_NB_NODES_PER_ELEMENT(_segment_2 ); break; } - case _segment_3 : { GET_NB_NODES_PER_ELEMENT(_segment_3 ); break; } - case _triangle_3 : { GET_NB_NODES_PER_ELEMENT(_triangle_3 ); break; } - case _triangle_6 : { GET_NB_NODES_PER_ELEMENT(_triangle_6 ); break; } - case _tetrahedron_4 : { GET_NB_NODES_PER_ELEMENT(_tetrahedron_4 ); break; } - case _tetrahedron_10 : { GET_NB_NODES_PER_ELEMENT(_tetrahedron_10); break; } - case _quadrangle_4 : { GET_NB_NODES_PER_ELEMENT(_quadrangle_4 ); break; } - case _point : { GET_NB_NODES_PER_ELEMENT(_point ); break; } - case _not_defined: - case _max_element_type: { - AKANTU_DEBUG_ERROR("Wrong type : " << type); - break; } - } - + AKANTU_BOOST_ELEMENT_SWITCH(GET_NB_NODES_PER_ELEMENT) #undef GET_NB_NODES_PER_ELEMENT AKANTU_DEBUG_OUT(); return nb_nodes_per_element; } /* -------------------------------------------------------------------------- */ inline ElementType Mesh::getP1ElementType(const ElementType & type) { AKANTU_DEBUG_IN(); ElementType element_p1 = _not_defined; #define GET_ELEMENT_P1(type) \ element_p1 = ElementClass::getP1ElementType() - switch(type) { - case _segment_2 : { GET_ELEMENT_P1(_segment_2 ); break; } - case _segment_3 : { GET_ELEMENT_P1(_segment_3 ); break; } - case _triangle_3 : { GET_ELEMENT_P1(_triangle_3 ); break; } - case _triangle_6 : { GET_ELEMENT_P1(_triangle_6 ); break; } - case _tetrahedron_4 : { GET_ELEMENT_P1(_tetrahedron_4 ); break; } - case _tetrahedron_10 : { GET_ELEMENT_P1(_tetrahedron_10); break; } - case _quadrangle_4 : { GET_ELEMENT_P1(_quadrangle_4 ); break; } - case _point : { GET_ELEMENT_P1(_point ); break; } - case _not_defined: - case _max_element_type: { - AKANTU_DEBUG_ERROR("Wrong type : " << type); - break; } - } - + AKANTU_BOOST_ELEMENT_SWITCH(GET_ELEMENT_P1) #undef GET_NB_NODES_PER_ELEMENT_P1 AKANTU_DEBUG_OUT(); return element_p1; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getSpatialDimension(const ElementType & type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = 0; #define GET_SPATIAL_DIMENSION(type) \ spatial_dimension = ElementClass::getSpatialDimension() - switch(type) { - case _segment_2 : { GET_SPATIAL_DIMENSION(_segment_2 ); break; } - case _segment_3 : { GET_SPATIAL_DIMENSION(_segment_3 ); break; } - case _triangle_3 : { GET_SPATIAL_DIMENSION(_triangle_3 ); break; } - case _triangle_6 : { GET_SPATIAL_DIMENSION(_triangle_6 ); break; } - case _tetrahedron_4 : { GET_SPATIAL_DIMENSION(_tetrahedron_4 ); break; } - case _tetrahedron_10 : { GET_SPATIAL_DIMENSION(_tetrahedron_10); break; } - case _quadrangle_4 : { GET_SPATIAL_DIMENSION(_quadrangle_4 ); break; } - case _point : { GET_SPATIAL_DIMENSION(_point ); break; } - case _not_defined: - case _max_element_type: { - AKANTU_DEBUG_ERROR("Wrong type : " << type); - break; } - } - + AKANTU_BOOST_ELEMENT_SWITCH(GET_SPATIAL_DIMENSION) #undef GET_SPATIAL_DIMENSION AKANTU_DEBUG_OUT(); return spatial_dimension; } /* -------------------------------------------------------------------------- */ inline ElementType Mesh::getFacetElementType(const ElementType & type) { AKANTU_DEBUG_IN(); ElementType surface_type = _not_defined; #define GET_FACET_TYPE(type) \ surface_type = ElementClass::getFacetElementType() - switch(type) { - case _segment_2 : { GET_FACET_TYPE(_segment_2 ); break; } - case _segment_3 : { GET_FACET_TYPE(_segment_3 ); break; } - case _triangle_3 : { GET_FACET_TYPE(_triangle_3 ); break; } - case _triangle_6 : { GET_FACET_TYPE(_triangle_6 ); break; } - case _tetrahedron_4 : { GET_FACET_TYPE(_tetrahedron_4 ); break; } - case _tetrahedron_10 : { GET_FACET_TYPE(_tetrahedron_10); break; } - case _quadrangle_4 : { GET_FACET_TYPE(_quadrangle_4 ); break; } - case _point : { GET_FACET_TYPE(_point ); break; } - case _not_defined: - case _max_element_type: { - AKANTU_DEBUG_ERROR("Wrong type : " << type); - break; } - } - + AKANTU_BOOST_ELEMENT_SWITCH(GET_FACET_TYPE) #undef GET_FACET_TYPE AKANTU_DEBUG_OUT(); return surface_type; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbFacetsPerElement(const ElementType & type) { AKANTU_DEBUG_IN(); UInt n_facet = 0; #define GET_NB_FACET(type) \ n_facet = ElementClass::getNbFacetsPerElement() - switch(type) { - case _segment_2 : { GET_NB_FACET(_segment_2 ); break; } - case _segment_3 : { GET_NB_FACET(_segment_3 ); break; } - case _triangle_3 : { GET_NB_FACET(_triangle_3 ); break; } - case _triangle_6 : { GET_NB_FACET(_triangle_6 ); break; } - case _tetrahedron_4 : { GET_NB_FACET(_tetrahedron_4 ); break; } - case _tetrahedron_10 : { GET_NB_FACET(_tetrahedron_10); break; } - case _quadrangle_4 : { GET_NB_FACET(_quadrangle_4 ); break; } - case _point : { GET_NB_FACET(_point ); break; } - case _not_defined: - case _max_element_type: { - AKANTU_DEBUG_ERROR("Wrong type : " << type); - break; } - } - + AKANTU_BOOST_ELEMENT_SWITCH(GET_NB_FACET) #undef GET_NB_FACET AKANTU_DEBUG_OUT(); return n_facet; } /* -------------------------------------------------------------------------- */ inline UInt ** Mesh::getFacetLocalConnectivity(const ElementType & type) { AKANTU_DEBUG_IN(); UInt ** facet_conn = NULL; #define GET_FACET_CON(type) \ facet_conn = ElementClass::getFacetLocalConnectivityPerElement() - switch(type) { - case _segment_2 : { GET_FACET_CON(_segment_2 ); break; } - case _segment_3 : { GET_FACET_CON(_segment_3 ); break; } - case _triangle_3 : { GET_FACET_CON(_triangle_3 ); break; } - case _triangle_6 : { GET_FACET_CON(_triangle_6 ); break; } - case _tetrahedron_4 : { GET_FACET_CON(_tetrahedron_4 ); break; } - case _tetrahedron_10 : { GET_FACET_CON(_tetrahedron_10); break; } - case _quadrangle_4 : { GET_FACET_CON(_quadrangle_4 ); break; } - case _point : { GET_FACET_CON(_point ); break; } - case _not_defined: - case _max_element_type: { - AKANTU_DEBUG_ERROR("Wrong type : " << type); - break; } - } - + AKANTU_BOOST_ELEMENT_SWITCH(GET_FACET_CON) #undef GET_FACET_CON AKANTU_DEBUG_OUT(); return facet_conn; } diff --git a/mesh_utils/mesh_pbc.cc b/mesh_utils/mesh_pbc.cc new file mode 100644 index 000000000..7c15e5bfb --- /dev/null +++ b/mesh_utils/mesh_pbc.cc @@ -0,0 +1,259 @@ +/** + * @file mesh_pbc.cc + * @author anciaux + * @date Tue Feb 8 10:48:01 2011 + * + * @brief periodic boundary condition connectivity tweak + * + * @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 "mesh_utils.hh" +#include + +__BEGIN_AKANTU__ + +/* -------------------------------------------------------------------------- */ +/// class that sorts a set of nodes of same coordinates in 'dir' direction +class CoordinatesComparison { +public: + CoordinatesComparison (const UInt dimension, + const UInt dirx, const UInt diry, + Real * coords): + dim(dimension),dir_x(dirx),dir_y(diry),coordinates(coords){} + + bool operator() (UInt n1, UInt n2){ + Real p1_x = coordinates[dim*n1+dir_x]; + Real p2_x = coordinates[dim*n2+dir_x]; + Real p1_y = coordinates[dim*n1+dir_y]; + Real p2_y = coordinates[dim*n2+dir_y]; + Real diff_x = p1_x - p2_x; + if (fabs(diff_x) > Math::tolerance) + return diff_x > 0.0 ? false : true; + else { + Real diff_y = p1_y - p2_y; + return diff_y > 0 ? false : true; + } + } +private: + UInt dim; + UInt dir_x; + UInt dir_y; + Real * coordinates; +}; + +/* -------------------------------------------------------------------------- */ +void MeshUtils::tweakConnectivityForPBC(Mesh & mesh, + UInt flag_x,UInt flag_y,UInt flag_z){ + std::map pbc_pair; + mesh.computeBoundingBox(); + + if (flag_x) computePBCMap(mesh,0,pbc_pair); + if (flag_y) computePBCMap(mesh,1,pbc_pair); + if (flag_z) computePBCMap(mesh,2,pbc_pair); + + { + std::map::iterator it = pbc_pair.begin(); + std::map::iterator end = pbc_pair.end(); + + Real * coords = mesh.nodes->values; + UInt dim = mesh.getSpatialDimension(); + while(it != end){ + UInt i1 = (*it).first; + UInt i2 = (*it).second; + + AKANTU_DEBUG_INFO("pairing " << i1 << "(" + << coords[dim*i1] << "," << coords[dim*i1+1] << "," << coords[dim*i1+2] + << ") with" + << i2 << "(" + << coords[dim*i2] << "," << coords[dim*i2+1] << "," << coords[dim*i2+2] + << ")"); + ++it; + } + } + // now loop over the elements to change the connectivity of some elements + const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); + Mesh::ConnectivityTypeList::const_iterator it; + for(it = type_list.begin(); it != type_list.end(); ++it) { + UInt nb_elem = mesh.getNbElement(*it); + UInt nb_node_per_elem = mesh.getNbNodesPerElement(*it); + UInt * conn = mesh.getConnectivityPointer(*it)->values; + for (UInt index = 0; index < nb_elem*nb_node_per_elem; index++) { + if (pbc_pair.count(conn[index])){ + conn[index] = pbc_pair[conn[index]]; + } + } + } + +} + + +/* -------------------------------------------------------------------------- */ + +void MeshUtils::computePBCMap(const Mesh & mymesh,const UInt dir, + std::map & pbc_pair){ + std::vector selected_left; + std::vector selected_right; + + Real * coords = mymesh.nodes->values; + const UInt nb_nodes = mymesh.nodes->getSize(); + const UInt dim = mymesh.getSpatialDimension(); + AKANTU_DEBUG_INFO("min " << mymesh.xmin[dir]); + AKANTU_DEBUG_INFO("max " << mymesh.xmax[dir]); + + for (UInt i = 0; i < nb_nodes; ++i) { + AKANTU_DEBUG_TRACE("treating " << coords[dim*i+dir]); + if(Math::are_float_equal(coords[dim*i+dir],mymesh.xmin[dir])){ + AKANTU_DEBUG_TRACE("pushing node " << i << " on the left side"); + selected_left.push_back(i); + } + else if(Math::are_float_equal(coords[dim*i+dir],mymesh.xmax[dir])){ + selected_right.push_back(i); + AKANTU_DEBUG_TRACE("pushing node " << i << " on the right side"); + } + } + + AKANTU_DEBUG_INFO("found " << selected_left.size() << " and " << selected_right.size() + << " nodes at each boundary for direction " << dir); + + UInt dir_x,dir_y; + + if (dir == 0){ + dir_x = 1;dir_y = 2; + } + else if (dir == 1){ + dir_x = 0;dir_y = 2; + } + else if (dir == 2){ + dir_x = 0;dir_y = 1; + } + + CoordinatesComparison compare_nodes(dim,dir_x,dir_y,coords); + + std::sort(selected_left.begin(),selected_left.end(),compare_nodes); + std::sort(selected_right.begin(),selected_right.end(),compare_nodes); + + std::vector::iterator it_left = selected_left.begin(); + std::vector::iterator end_left = selected_left.end(); + + std::vector::iterator it_right = selected_right.begin(); + std::vector::iterator end_right = selected_right.end(); + + while ((it_left != end_left) && (it_right != end_right) ){ + UInt i1 = *it_left; + UInt i2 = *it_right; + + AKANTU_DEBUG_TRACE("do I pair? " << i1 << "(" + << coords[dim*i1] << "," << coords[dim*i1+1] << "," << coords[dim*i1+2] + << ") with" + << i2 << "(" + << coords[dim*i2] << "," << coords[dim*i2+1] << "," << coords[dim*i2+2] + << ") in direction " << dir); + + + UInt dx = coords[dim*i1+dir_x] - coords[dim*i2+dir_x]; + UInt dy = coords[dim*i1+dir_y] - coords[dim*i2+dir_y]; + + if (fabs(dx*dx+dy*dy) < Math::tolerance) + { + //then i match these pairs + if (pbc_pair.count(i2)){ + i2 = pbc_pair[i2]; + } + pbc_pair[i1] = i2; + + AKANTU_DEBUG_TRACE("pairing " << i1 << "(" + << coords[dim*i1] << "," << coords[dim*i1+1] << "," << coords[dim*i1+2] + << ") with" + << i2 << "(" + << coords[dim*i2] << "," << coords[dim*i2+1] << "," << coords[dim*i2+2] + << ") in direction " << dir); + ++it_left; + ++it_right; + } + else if (fabs(dy) < Math::tolerance && dx > 0) ++it_right; + else if (fabs(dy) < Math::tolerance && dx < 0) ++it_left; + else if (dy > 0) ++it_right; + else if (dy < 0) ++it_left; + else { + AKANTU_DEBUG_ERROR("this should not append"); + } + } + AKANTU_DEBUG_INFO("found " << pbc_pair.size() << " pairs for direction " << dir); + +} + +// // now put together the map in such a way that only one indirection is associating pbc nodes +// int cpt=0; +// is_pair_at_corner = malloc(sizeof(int)*(n_pbc_pair_x+n_pbc_pair_y+n_pbc_pair_z)); +// memset(is_pair_at_corner,0,sizeof(int)*(n_pbc_pair_x+n_pbc_pair_y+n_pbc_pair_z)); +// if (pbc_nodes_1) free(pbc_nodes_1); +// pbc_nodes_1 = malloc(sizeof(int)*(n_pbc_pair_x+n_pbc_pair_y+n_pbc_pair_z)); +// if (pbc_nodes_2) free(pbc_nodes_2); +// pbc_nodes_2 = malloc(sizeof(int)*(n_pbc_pair_x+n_pbc_pair_y+n_pbc_pair_z)); + +// for(i=0; i + * @author Guillaume ANCIAUX * @date Wed Aug 18 14:21:00 2010 * * @brief All mesh utils necessary for various tasks * * @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 "mesh_utils.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ void MeshUtils::buildNode2Elements(const Mesh & mesh, Vector & node_offset, Vector & node_to_elem, UInt spatial_dimension) { AKANTU_DEBUG_IN(); if (spatial_dimension == 0) spatial_dimension = mesh.getSpatialDimension(); UInt nb_nodes = mesh.getNbNodes(); const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; UInt nb_types = type_list.size(); UInt nb_good_types = 0; UInt nb_nodes_per_element[nb_types]; UInt nb_nodes_per_element_p1[nb_types]; UInt * conn_val[nb_types]; UInt nb_element[nb_types]; ElementType type_p1; for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; if(Mesh::getSpatialDimension(type) != spatial_dimension) continue; nb_nodes_per_element[nb_good_types] = Mesh::getNbNodesPerElement(type); type_p1 = Mesh::getP1ElementType(type); nb_nodes_per_element_p1[nb_good_types] = Mesh::getNbNodesPerElement(type_p1); conn_val[nb_good_types] = mesh.getConnectivity(type).values; nb_element[nb_good_types] = mesh.getConnectivity(type).getSize(); nb_good_types++; } AKANTU_DEBUG_ASSERT(nb_good_types != 0, "Some elements must be found in right dimension to compute facets!"); /// array for the node-element list node_offset.resize(nb_nodes + 1); UInt * node_offset_val = node_offset.values; /// count number of occurrence of each node memset(node_offset_val, 0, (nb_nodes + 1)*sizeof(UInt)); for (UInt t = 0; t < nb_good_types; ++t) { for (UInt el = 0; el < nb_element[t]; ++el) { UInt el_offset = el*nb_nodes_per_element[t]; for (UInt n = 0; n < nb_nodes_per_element_p1[t]; ++n) { node_offset_val[conn_val[t][el_offset + n]]++; } } } /// convert the occurrence array in a csr one for (UInt i = 1; i < nb_nodes; ++i) node_offset_val[i] += node_offset_val[i-1]; for (UInt i = nb_nodes; i > 0; --i) node_offset_val[i] = node_offset_val[i-1]; node_offset_val[0] = 0; /// rearrange element to get the node-element list node_to_elem.resize(node_offset_val[nb_nodes]); UInt * node_to_elem_val = node_to_elem.values; for (UInt t = 0, linearized_el = 0; t < nb_good_types; ++t) for (UInt el = 0; el < nb_element[t]; ++el, ++linearized_el) { UInt el_offset = el*nb_nodes_per_element[t]; for (UInt n = 0; n < nb_nodes_per_element_p1[t]; ++n) node_to_elem_val[node_offset_val[conn_val[t][el_offset + n]]++] = linearized_el; } for (UInt i = nb_nodes; i > 0; --i) node_offset_val[i] = node_offset_val[i-1]; node_offset_val[0] = 0; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildNode2ElementsByElementType(const Mesh & mesh, ElementType type, Vector & node_offset, Vector & node_to_elem) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_elements = mesh.getConnectivity(type).getSize(); UInt * conn_val = mesh.getConnectivity(type).values; /// array for the node-element list node_offset.resize(nb_nodes + 1); UInt * node_offset_val = node_offset.values; memset(node_offset_val, 0, (nb_nodes + 1)*sizeof(UInt)); /// count number of occurrence of each node for (UInt el = 0; el < nb_elements; ++el) for (UInt n = 0; n < nb_nodes_per_element; ++n) node_offset_val[conn_val[nb_nodes_per_element*el + n]]++; /// convert the occurrence array in a csr one for (UInt i = 1; i < nb_nodes; ++i) node_offset_val[i] += node_offset_val[i-1]; for (UInt i = nb_nodes; i > 0; --i) node_offset_val[i] = node_offset_val[i-1]; node_offset_val[0] = 0; /// save the element index in the node-element list node_to_elem.resize(node_offset_val[nb_nodes]); UInt * node_to_elem_val = node_to_elem.values; for (UInt el = 0; el < nb_elements; ++el) for (UInt n = 0; n < nb_nodes_per_element; ++n) { node_to_elem_val[node_offset_val[conn_val[nb_nodes_per_element*el + n]]] = el; node_offset_val[conn_val[nb_nodes_per_element*el + n]]++; } /// rearrange node_offset to start with 0 for (UInt i = nb_nodes; i > 0; --i) node_offset_val[i] = node_offset_val[i-1]; node_offset_val[0] = 0; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildFacets(Mesh & mesh, bool boundary_flag, bool internal_flag){ AKANTU_DEBUG_IN(); Vector node_offset; Vector node_to_elem; buildNode2Elements(mesh,node_offset,node_to_elem); // std::cout << "node offset " << std::endl << node_offset << std::endl; // std::cout << "node to elem " << std::endl << node_to_elem << std::endl; const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; UInt nb_types = type_list.size(); UInt nb_good_types = 0; UInt nb_nodes_per_element[nb_types]; UInt nb_nodes_per_facet[nb_types]; UInt nb_facets[nb_types]; UInt ** node_in_facet[nb_types]; Vector * connectivity_facets[nb_types]; Vector * connectivity_internal_facets[nb_types]; UInt * conn_val[nb_types]; UInt nb_element[nb_types]; ElementType facet_type; Mesh * internal_facets_mesh = NULL; UInt spatial_dimension = mesh.getSpatialDimension(); if (internal_flag) internal_facets_mesh = mesh.getInternalFacetsMeshPointer(); for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; if(mesh.getSpatialDimension(type) != spatial_dimension) continue; nb_nodes_per_element[nb_good_types] = mesh.getNbNodesPerElement(type); facet_type = mesh.getFacetElementType(type); nb_facets[nb_good_types] = mesh.getNbFacetsPerElement(type); node_in_facet[nb_good_types] = mesh.getFacetLocalConnectivity(type); nb_nodes_per_facet[nb_good_types] = mesh.getNbNodesPerElement(facet_type); if (boundary_flag){ // getting connectivity of boundary facets connectivity_facets[nb_good_types] = mesh.getConnectivityPointer(facet_type); connectivity_facets[nb_good_types]->resize(0); } if (internal_flag){ // getting connectivity of internal facets connectivity_internal_facets[nb_good_types] = internal_facets_mesh->getConnectivityPointer(facet_type); connectivity_internal_facets[nb_good_types]->resize(0); } conn_val[nb_good_types] = mesh.getConnectivity(type).values; nb_element[nb_good_types] = mesh.getConnectivity(type).getSize(); nb_good_types++; } Vector counter; /// count number of occurrence of each node for (UInt t = 0,linearized_el = 0; t < nb_good_types; ++t) { for (UInt el = 0; el < nb_element[t]; ++el, ++linearized_el) { UInt el_offset = el*nb_nodes_per_element[t]; for (UInt f = 0; f < nb_facets[t]; ++f) { //build the nodes involved in facet 'f' UInt facet_nodes[nb_nodes_per_facet[t]]; for (UInt n = 0; n < nb_nodes_per_facet[t]; ++n) { UInt node_facet = node_in_facet[t][f][n]; facet_nodes[n] = conn_val[t][el_offset + node_facet]; } //our reference is the first node UInt * first_node_elements = node_to_elem.values+node_offset.values[facet_nodes[0]]; UInt first_node_nelements = node_offset.values[facet_nodes[0]+1]- node_offset.values[facet_nodes[0]]; counter.resize(first_node_nelements); memset(counter.values,0,sizeof(UInt)*first_node_nelements); //loop over the other nodes to search intersecting elements for (UInt n = 1; n < nb_nodes_per_facet[t]; ++n) { UInt * node_elements = node_to_elem.values+node_offset.values[facet_nodes[n]]; UInt node_nelements = node_offset.values[facet_nodes[n]+1]- node_offset.values[facet_nodes[n]]; for (UInt el1 = 0; el1 < first_node_nelements; ++el1) { for (UInt el2 = 0; el2 < node_nelements; ++el2) { if (first_node_elements[el1] == node_elements[el2]) { ++counter.values[el1]; break; } } if (counter.values[el1] == nb_nodes_per_facet[t]) break; } } bool connected_facet = false; //the connected elements are those for which counter is n_facet // UInt connected_element = -1; for (UInt el1 = 0; el1 < counter.getSize(); el1++) { UInt el_index = node_to_elem.values[node_offset.values[facet_nodes[0]]+el1]; if (counter.values[el1] == nb_nodes_per_facet[t]-1 && el_index > linearized_el){ // connected_element = el_index; AKANTU_DEBUG(dblDump,"connecting elements " << linearized_el << " and " << el_index); if (internal_flag) connectivity_internal_facets[t]->push_back(facet_nodes); } if (counter.values[el1] == nb_nodes_per_facet[t]-1 && el_index != linearized_el) connected_facet = true; } if (!connected_facet) { AKANTU_DEBUG(dblDump,"facet " << f << " in element " << linearized_el << " is a boundary"); if (boundary_flag) connectivity_facets[t]->push_back(facet_nodes); } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ // void MeshUtils::buildNormals(Mesh & mesh,UInt spatial_dimension){ // AKANTU_DEBUG_IN(); // const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); // Mesh::ConnectivityTypeList::const_iterator it; // UInt nb_types = type_list.size(); // UInt nb_nodes_per_element[nb_types]; // UInt nb_nodes_per_element_p1[nb_types]; // UInt nb_good_types = 0; // Vector * connectivity[nb_types]; // Vector * normals[nb_types]; // if (spatial_dimension == 0) spatial_dimension = mesh.getSpatialDimension(); // for(it = type_list.begin(); it != type_list.end(); ++it) { // ElementType type = *it; // ElementType type_p1 = mesh.getP1ElementType(type); // if(mesh.getSpatialDimension(type) != spatial_dimension) continue; // nb_nodes_per_element[nb_good_types] = mesh.getNbNodesPerElement(type); // nb_nodes_per_element_p1[nb_good_types] = mesh.getNbNodesPerElement(type_p1); // // getting connectivity // connectivity[nb_good_types] = mesh.getConnectivityPointer(type); // if (!connectivity[nb_good_types]) // AKANTU_DEBUG_ERROR("connectivity is not allocatted : this should probably not have happened"); // //getting array of normals // normals[nb_good_types] = mesh.getNormalsPointer(type); // if(normals[nb_good_types]) // normals[nb_good_types]->resize(0); // else // normals[nb_good_types] = &mesh.createNormals(type); // nb_good_types++; // } // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ void MeshUtils::renumberMeshNodes(Mesh & mesh, UInt * local_connectivities, UInt nb_local_element, UInt nb_ghost_element, ElementType type, Vector & nodes_numbers) { AKANTU_DEBUG_IN(); nodes_numbers.resize(0); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt * conn = local_connectivities; /// renumber the nodes for (UInt el = 0; el < nb_local_element + nb_ghost_element; ++el) { for (UInt n = 0; n < nb_nodes_per_element; ++n) { Int nn = nodes_numbers.find(*conn); if(nn == -1) { nodes_numbers.push_back(*conn); *conn = nodes_numbers.getSize() - 1; } else { *conn = nn; } *conn++; } } /// copy the renumbered connectivity to the right place Vector * local_conn = mesh.getConnectivityPointer(type); local_conn->resize(nb_local_element); memcpy(local_conn->values, local_connectivities, nb_local_element * nb_nodes_per_element * sizeof(UInt)); Vector * ghost_conn = mesh.getGhostConnectivityPointer(type); ghost_conn->resize(nb_ghost_element); memcpy(ghost_conn->values, local_connectivities + nb_local_element * nb_nodes_per_element, nb_ghost_element * nb_nodes_per_element * sizeof(UInt)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildSurfaceID(Mesh & mesh) { AKANTU_DEBUG_IN(); Vector node_offset; Vector node_to_elem; /// Get list of surface elements UInt spatial_dimension = mesh.getSpatialDimension(); buildNode2Elements(mesh, node_offset, node_to_elem, spatial_dimension-1); /// Find which types of elements have been linearized const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; UInt nb_types = type_list.size(); ElementType lin_element_type[nb_types]; UInt nb_lin_types = 0; UInt nb_nodes_per_element[nb_types]; UInt nb_nodes_per_element_p1[nb_types]; UInt * conn_val[nb_types]; UInt nb_element[nb_types]; ElementType type_p1; for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; if(Mesh::getSpatialDimension(type) != spatial_dimension) continue; ElementType facet_type = mesh.getFacetElementType(type); lin_element_type[nb_lin_types] = facet_type; nb_nodes_per_element[nb_lin_types] = Mesh::getNbNodesPerElement(facet_type); type_p1 = Mesh::getP1ElementType(facet_type); nb_nodes_per_element_p1[nb_lin_types] = Mesh::getNbNodesPerElement(type_p1); conn_val[nb_lin_types] = mesh.getConnectivity(facet_type).values; nb_element[nb_lin_types] = mesh.getConnectivity(facet_type).getSize(); nb_lin_types++; } for (UInt i = 1; i < nb_lin_types; ++i) nb_element[i] += nb_element[i+1]; for (UInt i = nb_lin_types; i > 0; --i) nb_element[i] = nb_element[i-1]; nb_element[0] = 0; /// Find close surfaces Vector surface_value_id(1, nb_element[nb_lin_types], -1); Int * surf_val = surface_value_id.values; UInt nb_surfaces = 0; UInt nb_cecked_elements; UInt nb_elements_to_ceck; UInt * elements_to_ceck = new UInt [nb_element[nb_lin_types]]; memset(elements_to_ceck, 0, nb_element[nb_lin_types]*sizeof(UInt)); for (UInt lin_el = 0; lin_el < nb_element[nb_lin_types]; ++lin_el) { if(surf_val[lin_el] != -1) continue; /* Surface id already assigned */ /* First element of new surface */ surf_val[lin_el] = nb_surfaces; nb_cecked_elements = 0; nb_elements_to_ceck = 1; memset(elements_to_ceck, 0, nb_element[nb_lin_types]*sizeof(UInt)); elements_to_ceck[0] = lin_el; // Find others elements belonging to this surface while(nb_cecked_elements < nb_elements_to_ceck) { UInt ceck_lin_el = elements_to_ceck[nb_cecked_elements]; // Transform linearized index of element into ElementType one UInt lin_type_nb = 0; while (ceck_lin_el >= nb_element[lin_type_nb+1]) lin_type_nb++; UInt ceck_el = ceck_lin_el - nb_element[lin_type_nb]; // Get connected elements UInt el_offset = ceck_el*nb_nodes_per_element[lin_type_nb]; for (UInt n = 0; n < nb_nodes_per_element_p1[lin_type_nb]; ++n) { UInt node_id = conn_val[lin_type_nb][el_offset + n]; for (UInt i = node_offset.values[node_id]; i < node_offset.values[node_id+1]; ++i) if(surf_val[node_to_elem.values[i]] == -1) { /* Found new surface element */ surf_val[node_to_elem.values[i]] = nb_surfaces; elements_to_ceck[nb_elements_to_ceck] = node_to_elem.values[i]; nb_elements_to_ceck++; } } nb_cecked_elements++; } nb_surfaces++; } delete [] elements_to_ceck; /// Transform local linearized element index in the global one for (UInt i = 0; i < nb_lin_types; ++i) nb_element[i] = nb_element[i+1] - nb_element[i]; UInt el_offset = 0; for (UInt type_it = 0; type_it < nb_lin_types; ++type_it) { ElementType type = lin_element_type[type_it]; Vector * surf_id_type = mesh.getSurfaceIdPointer(type); surf_id_type->resize(nb_element[type_it]); for (UInt el = 0; el < nb_element[type_it]; ++el) surf_id_type->values[el] = surf_val[el+el_offset]; el_offset += nb_element[type_it]; } /// Set nb_surfaces in mesh mesh.nb_surfaces = nb_surfaces; AKANTU_DEBUG_OUT(); } __END_AKANTU__ // LocalWords: ElementType diff --git a/mesh_utils/mesh_utils.hh b/mesh_utils/mesh_utils.hh index 99902bcab..208be7e74 100644 --- a/mesh_utils/mesh_utils.hh +++ b/mesh_utils/mesh_utils.hh @@ -1,105 +1,113 @@ /** * @file mesh_utils.hh - * @author Guillaume ANCIAUX + * @author Guillaume ANCIAUX * @date Wed Aug 18 14:03:39 2010 * * @brief All mesh utils necessary for various tasks * * @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_MESH_UTILS_HH__ #define __AKANTU_MESH_UTILS_HH__ #include "aka_common.hh" #include "mesh.hh" __BEGIN_AKANTU__ class MeshUtils { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MeshUtils(); virtual ~MeshUtils(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// build map from nodes to elements static void buildNode2Elements(const Mesh & mesh, Vector & node_offset, Vector & node_to_elem, UInt spatial_dimension = 0); /// build map from nodes to elements for a specific element type static void buildNode2ElementsByElementType(const Mesh & mesh, ElementType type, Vector & node_offset, Vector & node_to_elem); /// build facets elements : boundary and/or internals static void buildFacets(Mesh & mesh, bool boundary_flag=1, bool internal_flag=0); /// build normal to some elements // static void buildNormals(Mesh & mesh, UInt spatial_dimension=0); /// take the local_connectivity array as the array of local and ghost /// connectivity, renumber the nodes and set the connectivity of the mesh static void renumberMeshNodes(Mesh & mesh, UInt * local_connectivities, UInt nb_local_element, UInt nb_ghost_element, ElementType type, Vector & old_nodes); /// Detect closed surfaces of the mesh and save the surface id /// of the surface elements in the array surface_id static void buildSurfaceID(Mesh & mesh); + /// tweak mesh connectivity to activate pbc + static void tweakConnectivityForPBC(Mesh & mesh,UInt flag_x, + UInt flag_y=false,UInt flag_z=false); + /// function to print the contain of the class // virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: + /// compute pbc pair for on egiven direction + static void computePBCMap(const Mesh & mymesh,const UInt dir, + std::map & pbc_pair); + }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "mesh_utils_inline_impl.cc" /// standard output stream operator // inline std::ostream & operator <<(std::ostream & stream, const MeshUtils & _this) // { // _this.printself(stream); // return stream; // } __END_AKANTU__ #endif /* __AKANTU_MESH_UTILS_HH__ */ diff --git a/mesh_utils/mesh_utils_inline_impl.cc b/mesh_utils/mesh_utils_inline_impl.cc index 8e22eb91e..18b590f9b 100644 --- a/mesh_utils/mesh_utils_inline_impl.cc +++ b/mesh_utils/mesh_utils_inline_impl.cc @@ -1,30 +1,30 @@ /** * @file mesh_utils_inline_impl.cc - * @author Guillaume ANCIAUX + * @author Guillaume ANCIAUX * @date Wed Aug 18 14:05:36 2010 * * @brief All mesh utils necessary for various tasks * * @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 . * */ /* -------------------------------------------------------------------------- */ diff --git a/model/model_inline_impl.cc b/model/model_inline_impl.cc index 8b7495ad2..151f2609c 100644 --- a/model/model_inline_impl.cc +++ b/model/model_inline_impl.cc @@ -1,64 +1,64 @@ /** * @file model_inline_impl.cc - * @author Guillaume ANCIAUX + * @author Guillaume ANCIAUX * @date Fri Aug 20 17:18:08 2010 * * @brief inline implementation of the model 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 . * */ /* -------------------------------------------------------------------------- */ inline Model::Model(Mesh & mesh, UInt spatial_dimension, const ModelID & id, const MemoryID & memory_id) : Memory(memory_id), id(id) { AKANTU_DEBUG_IN(); this->spatial_dimension = (spatial_dimension == 0) ? mesh.getSpatialDimension() : spatial_dimension; std::stringstream sstr; sstr << id << ":fem"; this->fem = new FEM(mesh, spatial_dimension, sstr.str(), memory_id); this->fem_boundary = NULL; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline Model::~Model() { AKANTU_DEBUG_IN(); delete fem; if(fem_boundary) delete fem_boundary; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline FEM & Model::getFEMBoundary(){ AKANTU_DEBUG_IN(); if (!fem_boundary){ MeshUtils::buildFacets(fem->getMesh()); std::stringstream sstr; sstr << id << ":femboundary"; this->fem_boundary = new FEM(fem->getMesh(), spatial_dimension-1, sstr.str(), memory_id); } AKANTU_DEBUG_OUT(); return *this->fem_boundary; } diff --git a/model/solid_mechanics/material_parser.hh b/model/solid_mechanics/material_parser.hh index ff0dcde8a..66e9212d0 100644 --- a/model/solid_mechanics/material_parser.hh +++ b/model/solid_mechanics/material_parser.hh @@ -1,85 +1,85 @@ /** * @file material_parser.hh - * @author Guillaume ANCIAUX + * @author Guillaume ANCIAUX * @date Thu Nov 25 11:43:48 2010 * * @brief * * @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_MATERIAL_PARSER_HH__ #define __AKANTU_MATERIAL_PARSER_HH__ __BEGIN_AKANTU__ class MaterialParser { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialParser() : current_line(0) {}; virtual ~MaterialParser(){ infile.close(); }; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// open a file to parse void open(const std::string & filename); /// read the file and return the next material type std::string getNextMaterialType(); /// read properties and instanciate a given material object template Material * readMaterialObject(SolidMechanicsModel & model, MaterialID & mat_id); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: inline void my_getline(); /// number of the current line UInt current_line; /// material file std::ifstream infile; /// current read line std::string line; }; #include "material_parser_inline_impl.cc" __END_AKANTU__ #endif diff --git a/model/solid_mechanics/solid_mechanics_model_boundary.cc b/model/solid_mechanics/solid_mechanics_model_boundary.cc index a87fc9111..2912a4db6 100644 --- a/model/solid_mechanics/solid_mechanics_model_boundary.cc +++ b/model/solid_mechanics/solid_mechanics_model_boundary.cc @@ -1,197 +1,197 @@ /** * @file solid_mechanics_model_boundary.cc - * @author Guillaume ANCIAUX + * @author Guillaume ANCIAUX * @date Fri Nov 19 10:23:03 2010 * * @brief * * @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 "solid_mechanics_model.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /** * @param myf pointer to a function that fills a vector/tensor with respect to * passed coordinates * @param function_type flag to specify the take of function: _bft_stress is * tensor like and _bft_forces is traction like. */ void SolidMechanicsModel::computeForcesFromFunction(BoundaryFunction myf, BoundaryFunctionType function_type){ /** function type is ** _bft_forces : traction function is given ** _bft_stress : stress function is given */ std::stringstream name; name << id << ":solidmechanics:imposed_traction"; Vector traction_funct(0, spatial_dimension,name.str()); name.clear(); name << id << ":solidmechanics:imposed_stresses"; Vector stress_funct(0, spatial_dimension*spatial_dimension,name.str()); name.clear(); name << id << ":solidmechanics:quad_coords"; Vector quad_coords(0,spatial_dimension,"quad_coords"); UInt offset = 0; switch(function_type) { case _bft_stress: offset = spatial_dimension * spatial_dimension; break; case _bft_forces: offset = spatial_dimension; break; } //prepare the loop over element types const Mesh::ConnectivityTypeList & type_list = fem_boundary->getMesh().getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != fem_boundary->getElementDimension()) continue; UInt nb_quad = FEM::getNbQuadraturePoints(*it); UInt nb_element = fem_boundary->getMesh().getNbElement(*it); fem_boundary->interpolateOnQuadraturePoints(fem_boundary->getMesh().getNodes(), quad_coords, spatial_dimension, (*it)); Real * imposed_val = NULL; switch(function_type) { case _bft_stress: stress_funct.resize(nb_element*nb_quad); imposed_val = stress_funct.values; break; case _bft_forces: traction_funct.resize(nb_element*nb_quad); imposed_val = traction_funct.values; break; } /// sigma/tractions on each quadrature points Real * qcoord = quad_coords.values; for (UInt el = 0; el < nb_element; ++el) { for (UInt q = 0; q < nb_quad; ++q) { myf(qcoord, imposed_val); imposed_val += offset; qcoord += spatial_dimension; } } switch(function_type) { case _bft_stress: computeForcesByStressTensor(stress_funct,(*it)); break; case _bft_forces: computeForcesByTractionVector(traction_funct,(*it)); break; } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeForcesByStressTensor(const Vector & stresses, const ElementType & type){ AKANTU_DEBUG_IN(); UInt nb_element = fem_boundary->getMesh().getNbElement(type); UInt nb_quad = FEM::getNbQuadraturePoints(type); // check dimension match AKANTU_DEBUG_ASSERT(Mesh::getSpatialDimension(type) == fem_boundary->getElementDimension(), "element type dimension does not match the dimension of boundaries : " << fem_boundary->getElementDimension() << " != " << Mesh::getSpatialDimension(type)); // check size of the vector AKANTU_DEBUG_ASSERT(stresses.getSize() == nb_quad*nb_element, "the size of the vector should be the total number of quadrature points"); // check number of components AKANTU_DEBUG_ASSERT(stresses.getNbComponent() == spatial_dimension*spatial_dimension, "the number of components should be the dimension of 2-tensors"); std::stringstream name; name << id << ":solidmechanics:" << type << ":traction_boundary"; Vector funct(nb_element*nb_quad, spatial_dimension,name.str()); const Vector & normals_on_quad = fem_boundary->getNormalsOnQuadPoints(type); Math::matrix_vector(spatial_dimension,spatial_dimension,stresses,normals_on_quad,funct); computeForcesByTractionVector(funct,type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeForcesByTractionVector(const Vector & tractions, const ElementType & type){ AKANTU_DEBUG_IN(); UInt nb_element = fem_boundary->getMesh().getNbElement(type); UInt nb_nodes_per_element = fem_boundary->getMesh().getNbNodesPerElement(type); UInt nb_quad = FEM::getNbQuadraturePoints(type); // check dimension match AKANTU_DEBUG_ASSERT(Mesh::getSpatialDimension(type) == fem_boundary->getElementDimension(), "element type dimension does not match the dimension of boundaries : " << fem_boundary->getElementDimension() << " != " << Mesh::getSpatialDimension(type)); // check size of the vector AKANTU_DEBUG_ASSERT(tractions.getSize() == nb_quad*nb_element, "the size of the vector should be the total number of quadrature points"); // check number of components AKANTU_DEBUG_ASSERT(tractions.getNbComponent() == spatial_dimension, "the number of components should be the spatial dimension of the problem"); // do a complete copy of the vector Vector funct(tractions,true); // extend the vector to multiply by the shapes (prepare to assembly) funct.extendComponentsInterlaced(nb_nodes_per_element,spatial_dimension); // multiply by the shapes Real * funct_val = funct.values; Real * shapes_val = (fem_boundary->getShapes(type)).values; for (UInt el = 0; el < nb_element; ++el) { for (UInt q = 0; q < nb_quad; ++q) { for (UInt n = 0; n < nb_nodes_per_element; ++n,++shapes_val) { for (UInt i = 0; i < spatial_dimension; ++i) { *funct_val++ *= *shapes_val; } } } } // allocate the vector that will contain the integrated values std::stringstream name; name << id << ":solidmechanics:" << type << ":integral_boundary"; Vector int_funct(nb_element, spatial_dimension*nb_nodes_per_element,name.str()); //do the integration fem_boundary->integrate(funct, int_funct, spatial_dimension*nb_nodes_per_element, type); // assemble the result into force vector fem_boundary->assembleVector(int_funct,const_cast &>(getForce()), spatial_dimension, type); AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/model/solid_mechanics/solid_mechanics_model_material.cc b/model/solid_mechanics/solid_mechanics_model_material.cc index 320bd2985..1a8f64d8d 100644 --- a/model/solid_mechanics/solid_mechanics_model_material.cc +++ b/model/solid_mechanics/solid_mechanics_model_material.cc @@ -1,99 +1,99 @@ /** * @file solid_mechanics_model_material.cc - * @author Guillaume ANCIAUX + * @author Guillaume ANCIAUX * @date Thu Nov 25 10:48:53 2010 * * @brief * * @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 "solid_mechanics_model.hh" #include "material.hh" #include "aka_math.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::readMaterials(const std::string & filename) { MaterialParser parser; parser.open(filename.c_str()); std::string mat_type = parser.getNextMaterialType(); UInt mat_count = 0; while (mat_type != ""){ std::stringstream sstr_mat; sstr_mat << id << ":" << mat_count++ << ":" << mat_type; Material * material; MaterialID mat_id = sstr_mat.str(); /// read the material properties if(mat_type == "elastic") material = parser.readMaterialObject(*this,mat_id); else AKANTU_DEBUG_ERROR("Malformed material file : unknown material type " << mat_type); materials.push_back(material); mat_type = parser.getNextMaterialType(); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initMaterials() { AKANTU_DEBUG_ASSERT(materials.size() != 0, "No material to initialize !"); Material ** mat_val = &(materials.at(0)); /// fill the element filters of the materials using the element_material arrays const Mesh::ConnectivityTypeList & type_list = fem->getMesh().getConnectivityTypeList(_not_ghost); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; UInt nb_element = fem->getMesh().getNbElement(*it); UInt * elem_mat_val = element_material[*it]->values; for (UInt el = 0; el < nb_element; ++el) { mat_val[elem_mat_val[el]]->addElement(*it, el); } } /// @todo synchronize element material /// fill the element filters of the materials using the element_material arrays const Mesh::ConnectivityTypeList & ghost_type_list = fem->getMesh().getConnectivityTypeList(_ghost); for(it = ghost_type_list.begin(); it != ghost_type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; UInt nb_element = fem->getMesh().getNbGhostElement(*it); UInt * elem_mat_val = ghost_element_material[*it]->values; for (UInt el = 0; el < nb_element; ++el) { mat_val[elem_mat_val[el]]->addGhostElement(*it, el); } } std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { /// init internals properties (*mat_it)->initMaterial(); } } __END_AKANTU__ diff --git a/test/test_mesh_utils/CMakeLists.txt b/test/test_mesh_utils/CMakeLists.txt index 6438fa814..cddae74e9 100644 --- a/test/test_mesh_utils/CMakeLists.txt +++ b/test/test_mesh_utils/CMakeLists.txt @@ -1,38 +1,39 @@ #=============================================================================== # @file CMakeLists.txt # @author Nicolas Richart # @date Wed Feb 2 14:46:59 2011 # # @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 # #=============================================================================== #=============================================================================== # List of tests #=============================================================================== add_akantu_test(test_mesh_io "Test mesh io object") add_akantu_test(test_facet_extraction "Test mesh utils facet extraction") +add_akantu_test(test_pbc_tweak "Test pbc facilities") if(AKANTU_SCOTCH_ON) add_akantu_test(test_mesh_partitionate "Test mesh partition creation") endif(AKANTU_SCOTCH_ON) diff --git a/test/test_mesh_utils/test_facet_extraction/test_facet_extraction_tetrahedron_4.cc b/test/test_mesh_utils/test_facet_extraction/test_facet_extraction_tetrahedron_4.cc index 9519c546a..49a5f29c3 100644 --- a/test/test_mesh_utils/test_facet_extraction/test_facet_extraction_tetrahedron_4.cc +++ b/test/test_mesh_utils/test_facet_extraction/test_facet_extraction_tetrahedron_4.cc @@ -1,87 +1,87 @@ /** * @file test_facet_extraction_tetra1.cc - * @author Guillaume ANCIAUX + * @author Guillaume ANCIAUX * @date Thu Aug 19 13:05:27 2010 * * @brief test of internal facet extraction * * @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 "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER using namespace akantu; int main(int argc, char *argv[]) { int dim = 3; Mesh mesh(dim); MeshIOMSH mesh_io; mesh_io.read("cube.msh", mesh); MeshUtils::buildFacets(mesh,1,1); #ifdef AKANTU_USE_IOHELPER unsigned int nb_nodes = mesh.getNbNodes(); DumperParaview dumper; dumper.SetMode(TEXT); dumper.SetPoints(mesh.getNodes().values, dim, nb_nodes, "test-facet-extraction"); dumper.SetConnectivity((int*)mesh.getConnectivity(_tetrahedron_4).values, TETRA1, mesh.getNbElement(_tetrahedron_4), C_MODE); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); DumperParaview dumper_facet; dumper_facet.SetMode(TEXT); dumper_facet.SetPoints(mesh.getNodes().values, dim, nb_nodes, "test-facet-extraction_boundary"); dumper_facet.SetConnectivity((int*)mesh.getConnectivity(_triangle_3).values, TRIANGLE1, mesh.getNbElement(_triangle_3), C_MODE); dumper_facet.SetPrefix("paraview/"); dumper_facet.Init(); dumper_facet.Dump(); dumper_facet.SetPoints(mesh.getNodes().values, dim, nb_nodes, "test-facet-extraction_internal"); dumper_facet.SetConnectivity((int*)mesh.getInternalFacetsMesh().getConnectivity(_triangle_3).values, TRIANGLE1, mesh.getInternalFacetsMesh().getNbElement(_triangle_3), C_MODE); dumper_facet.Init(); dumper_facet.Dump(); #endif //AKANTU_USE_IOHELPER return EXIT_SUCCESS; } diff --git a/test/test_mesh_utils/test_facet_extraction/test_facet_extraction_triangle_3.cc b/test/test_mesh_utils/test_facet_extraction/test_facet_extraction_triangle_3.cc index 45b67e78e..af6d750d9 100644 --- a/test/test_mesh_utils/test_facet_extraction/test_facet_extraction_triangle_3.cc +++ b/test/test_mesh_utils/test_facet_extraction/test_facet_extraction_triangle_3.cc @@ -1,86 +1,86 @@ /** * @file test_facet_extraction.cc - * @author Guillaume ANCIAUX + * @author Guillaume ANCIAUX * @date Thu Aug 19 13:05:27 2010 * * @brief test of internal facet extraction * * @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 "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER using namespace akantu; int main(int argc, char *argv[]) { int dim = 2; Mesh mesh(dim); MeshIOMSH mesh_io; mesh_io.read("square.msh", mesh); MeshUtils::buildFacets(mesh,1,1); #ifdef AKANTU_USE_IOHELPER unsigned int nb_nodes = mesh.getNbNodes(); DumperParaview dumper; dumper.SetMode(TEXT); dumper.SetPoints(mesh.getNodes().values, dim, nb_nodes, "test-facet-extraction"); dumper.SetConnectivity((int*)mesh.getConnectivity(_triangle_3).values, TRIANGLE1, mesh.getNbElement(_triangle_3), C_MODE); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); DumperParaview dumper_facet; dumper_facet.SetMode(TEXT); dumper_facet.SetPoints(mesh.getNodes().values, dim, nb_nodes, "test-facet-extraction_boundary"); dumper_facet.SetConnectivity((int*)mesh.getConnectivity(_segment_2).values, LINE1, mesh.getNbElement(_segment_2), C_MODE); dumper_facet.SetPrefix("paraview/"); dumper_facet.Init(); dumper_facet.Dump(); dumper_facet.SetPoints(mesh.getNodes().values, dim, nb_nodes, "test-facet-extraction_internal"); dumper_facet.SetConnectivity((int*)mesh.getInternalFacetsMesh().getConnectivity(_segment_2).values, LINE1, mesh.getInternalFacetsMesh().getNbElement(_segment_2), C_MODE); dumper_facet.Init(); dumper_facet.Dump(); #endif //AKANTU_USE_IOHELPER return EXIT_SUCCESS; } diff --git a/test/test_mesh_utils/CMakeLists.txt b/test/test_mesh_utils/test_pbc_tweak/CMakeLists.txt similarity index 72% copy from test/test_mesh_utils/CMakeLists.txt copy to test/test_mesh_utils/test_pbc_tweak/CMakeLists.txt index 6438fa814..5f2974e39 100644 --- a/test/test_mesh_utils/CMakeLists.txt +++ b/test/test_mesh_utils/test_pbc_tweak/CMakeLists.txt @@ -1,38 +1,34 @@ #=============================================================================== # @file CMakeLists.txt -# @author Nicolas Richart -# @date Wed Feb 2 14:46:59 2011 +# @author Guillaume Anciaux +# @date Fri Jun 11 09:46:59 2010 # # @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 # #=============================================================================== #=============================================================================== -# List of tests -#=============================================================================== -add_akantu_test(test_mesh_io "Test mesh io object") -add_akantu_test(test_facet_extraction "Test mesh utils facet extraction") +add_mesh(test_pbc_cube_mesh cube.geo 3 1) -if(AKANTU_SCOTCH_ON) - add_akantu_test(test_mesh_partitionate "Test mesh partition creation") -endif(AKANTU_SCOTCH_ON) +register_test(test_pbc_tweak test_pbc_tweak.cc) +add_dependencies(test_pbc_tweak test_pbc_cube_mesh) diff --git a/test/test_mesh_utils/test_pbc_tweak/cube.geo b/test/test_mesh_utils/test_pbc_tweak/cube.geo new file mode 100644 index 000000000..f9014f701 --- /dev/null +++ b/test/test_mesh_utils/test_pbc_tweak/cube.geo @@ -0,0 +1,132 @@ +h = 1; + +Point(1) = {0.0, 0.0, -1.0, h}; +Point(2) = {1.0, 0.0, -1.0, h}; +Point(3) = {0.0, 1.0, -1.0, h}; +Point(4) = {1.0, 1.0, -1.0, h}; +Point(5) = {-1.0,0.0, -1.0, h}; +Point(6) = {0.0,-1.0, -1.0, h}; +Point(7) = {-1.0,-1.0, -1.0, h}; +Point(8) = {1.0, -1.0, -1.0, h}; +Point(9) = {-1.0, 1.0, -1.0, h}; + +Point(11) = {0.0, 0.0, 0.0, h}; +Point(12) = {1.0, 0.0, 0.0, h}; +Point(13) = {0.0, 1.0, 0.0, h}; +Point(14) = {1.0, 1.0, 0.0, h}; +Point(15) = {-1.0,0.0, 0.0, h}; +Point(16) = {0.0,-1.0, 0.0, h}; +Point(17) = {-1.0,-1.0, 0.0, h}; +Point(18) = {1.0, -1.0, 0.0, h}; +Point(19) = {-1.0, 1.0, 0.0, h}; + +Point(21) = {0.0, 0.0, 1.0, h}; +Point(22) = {1.0, 0.0, 1.0, h}; +Point(23) = {0.0, 1.0, 1.0, h}; +Point(24) = {1.0, 1.0, 1.0, h}; +Point(25) = {-1.0,0.0, 1.0, h}; +Point(26) = {0.0,-1.0, 1.0, h}; +Point(27) = {-1.0,-1.0, 1.0, h}; +Point(28) = {1.0, -1.0, 1.0, h}; +Point(29) = {-1.0, 1.0, 1.0, h}; + +Line(1) = {13, 23}; +Line(2) = {23, 24}; +Line(3) = {24, 14}; +Line(4) = {14, 13}; +Line(5) = {23, 29}; +Line(6) = {29, 19}; +Line(7) = {19, 13}; +Line(8) = {14, 13}; +Line(9) = {13, 3}; +Line(10) = {3, 4}; +Line(11) = {4, 14}; +Line(12) = {19, 9}; +Line(13) = {9, 3}; + +Line(14) = {29, 25}; +Line(15) = {25, 21}; +Line(16) = {21, 23}; +Line(17) = {21, 22}; +Line(18) = {22, 24}; +Line(19) = {25, 27}; +Line(20) = {27, 26}; +Line(21) = {26, 21}; +Line(22) = {26, 28}; +Line(23) = {28, 22}; +Line(24) = {22, 12}; +Line(25) = {12, 14}; +Line(26) = {12, 2}; +Line(27) = {2, 4}; +Line(28) = {28, 18}; +Line(29) = {18, 12}; +Line(30) = {18, 8}; +Line(31) = {8, 2}; +Line(32) = {2, 1}; +Line(33) = {1, 3}; +Line(34) = {8, 6}; +Line(35) = {6, 1}; +Line(36) = {1, 5}; +Line(37) = {5, 9}; +Line(38) = {5, 7}; +Line(39) = {7, 6}; +Line(40) = {5, 15}; +Line(41) = {15, 19}; +Line(42) = {15, 25}; +Line(43) = {7, 17}; +Line(44) = {17, 15}; +Line(45) = {17, 27}; +Line(46) = {6, 16}; +Line(47) = {16, 17}; +Line(48) = {16, 26}; +Line(49) = {16, 18}; +Line Loop(50) = {5, 14, 15, 16}; +Plane Surface(51) = {50}; +Line Loop(52) = {2, -18, -17, 16}; +Plane Surface(53) = {52}; +Line Loop(54) = {15, -21, -20, -19}; +Plane Surface(55) = {54}; +Line Loop(56) = {17, -23, -22, 21}; +Plane Surface(57) = {56}; +Line Loop(58) = {3, -25, -24, 18}; +Plane Surface(59) = {58}; +Line Loop(60) = {11, -25, 26, 27}; +Plane Surface(61) = {60}; +Line Loop(62) = {24, -29, -28, 23}; +Plane Surface(63) = {62}; +Line Loop(64) = {26, -31, -30, 29}; +Plane Surface(65) = {64}; +Line Loop(66) = {33, 10, -27, 32}; +Plane Surface(67) = {66}; +Line Loop(68) = {32, -35, -34, 31}; +Plane Surface(69) = {68}; +Line Loop(70) = {36, 38, 39, 35}; +Plane Surface(71) = {70}; +Line Loop(72) = {37, 13, -33, 36}; +Plane Surface(73) = {72}; +Line Loop(74) = {14, -42, 41, -6}; +Plane Surface(75) = {74}; +Line Loop(76) = {19, -45, 44, 42}; +Plane Surface(77) = {76}; +Line Loop(78) = {41, 12, -37, 40}; +Plane Surface(79) = {78}; +Line Loop(80) = {44, -40, 38, 43}; +Plane Surface(81) = {80}; +Line Loop(82) = {2, 3, 4, 1}; +Plane Surface(83) = {82}; +Line Loop(84) = {4, 9, 10, 11}; +Plane Surface(85) = {84}; +Line Loop(86) = {5, 6, 7, 1}; +Plane Surface(87) = {86}; +Line Loop(88) = {7, 9, -13, -12}; +Plane Surface(89) = {88}; +Line Loop(90) = {20, -48, 47, 45}; +Plane Surface(91) = {90}; +Line Loop(92) = {47, -43, 39, 46}; +Plane Surface(93) = {92}; +Line Loop(94) = {22, 28, -49, 48}; +Plane Surface(95) = {94}; +Line Loop(96) = {49, 30, 34, 46}; +Plane Surface(97) = {96}; +Surface Loop(98) = {77, 55, 51, 87, 75, 79, 89, 85, 83, 53, 59, 61, 65, 69, 67, 73, 71, 81, 93, 91, 95, 57, 63, 97}; +Volume(99) = {98}; diff --git a/test/test_mesh_utils/test_facet_extraction/test_facet_extraction_tetrahedron_4.cc b/test/test_mesh_utils/test_pbc_tweak/test_pbc_tweak.cc similarity index 55% copy from test/test_mesh_utils/test_facet_extraction/test_facet_extraction_tetrahedron_4.cc copy to test/test_mesh_utils/test_pbc_tweak/test_pbc_tweak.cc index 9519c546a..51d8700a9 100644 --- a/test/test_mesh_utils/test_facet_extraction/test_facet_extraction_tetrahedron_4.cc +++ b/test/test_mesh_utils/test_pbc_tweak/test_pbc_tweak.cc @@ -1,87 +1,87 @@ /** * @file test_facet_extraction_tetra1.cc - * @author Guillaume ANCIAUX + * @author Guillaume ANCIAUX * @date Thu Aug 19 13:05:27 2010 * * @brief test of internal facet extraction * * @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 "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER using namespace akantu; int main(int argc, char *argv[]) { int dim = 3; Mesh mesh(dim); MeshIOMSH mesh_io; mesh_io.read("cube.msh", mesh); - MeshUtils::buildFacets(mesh,1,1); + MeshUtils::tweakConnectivityForPBC(mesh,1,1,1); -#ifdef AKANTU_USE_IOHELPER - unsigned int nb_nodes = mesh.getNbNodes(); - DumperParaview dumper; - dumper.SetMode(TEXT); +// #ifdef AKANTU_USE_IOHELPER +// unsigned int nb_nodes = mesh.getNbNodes(); +// DumperParaview dumper; +// dumper.SetMode(TEXT); - dumper.SetPoints(mesh.getNodes().values, dim, nb_nodes, "test-facet-extraction"); - dumper.SetConnectivity((int*)mesh.getConnectivity(_tetrahedron_4).values, - TETRA1, mesh.getNbElement(_tetrahedron_4), C_MODE); - dumper.SetPrefix("paraview/"); - dumper.Init(); - dumper.Dump(); +// dumper.SetPoints(mesh.getNodes().values, dim, nb_nodes, "test-facet-extraction"); +// dumper.SetConnectivity((int*)mesh.getConnectivity(_tetrahedron_4).values, +// TETRA1, mesh.getNbElement(_tetrahedron_4), C_MODE); +// dumper.SetPrefix("paraview/"); +// dumper.Init(); +// dumper.Dump(); - DumperParaview dumper_facet; - dumper_facet.SetMode(TEXT); +// DumperParaview dumper_facet; +// dumper_facet.SetMode(TEXT); - dumper_facet.SetPoints(mesh.getNodes().values, dim, nb_nodes, "test-facet-extraction_boundary"); - dumper_facet.SetConnectivity((int*)mesh.getConnectivity(_triangle_3).values, - TRIANGLE1, mesh.getNbElement(_triangle_3), C_MODE); - dumper_facet.SetPrefix("paraview/"); - dumper_facet.Init(); - dumper_facet.Dump(); +// dumper_facet.SetPoints(mesh.getNodes().values, dim, nb_nodes, "test-facet-extraction_boundary"); +// dumper_facet.SetConnectivity((int*)mesh.getConnectivity(_triangle_3).values, +// TRIANGLE1, mesh.getNbElement(_triangle_3), C_MODE); +// dumper_facet.SetPrefix("paraview/"); +// dumper_facet.Init(); +// dumper_facet.Dump(); - dumper_facet.SetPoints(mesh.getNodes().values, dim, nb_nodes, "test-facet-extraction_internal"); - dumper_facet.SetConnectivity((int*)mesh.getInternalFacetsMesh().getConnectivity(_triangle_3).values, - TRIANGLE1, mesh.getInternalFacetsMesh().getNbElement(_triangle_3), C_MODE); - dumper_facet.Init(); - dumper_facet.Dump(); +// dumper_facet.SetPoints(mesh.getNodes().values, dim, nb_nodes, "test-facet-extraction_internal"); +// dumper_facet.SetConnectivity((int*)mesh.getInternalFacetsMesh().getConnectivity(_triangle_3).values, +// TRIANGLE1, mesh.getInternalFacetsMesh().getNbElement(_triangle_3), C_MODE); +// dumper_facet.Init(); +// dumper_facet.Dump(); -#endif //AKANTU_USE_IOHELPER +//#endif //AKANTU_USE_IOHELPER return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_cube3d.cc b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_cube3d.cc index 55dbdc28e..a4fc2dfdd 100644 --- a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_cube3d.cc +++ b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_cube3d.cc @@ -1,140 +1,140 @@ /** * @file test_solid_mechanics_model_cube3d.cc - * @author Guillaume ANCIAUX + * @author Guillaume ANCIAUX * @date Tue Aug 17 11:31:22 2010 * * @brief test of the class SolidMechanicsModel on the 3d cube * * @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 "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "solid_mechanics_model.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER int main(int argc, char *argv[]) { akantu::UInt max_steps = 10000; akantu::Real epot, ekin; #ifdef AKANTU_USE_IOHELPER akantu::ElementType type = akantu::_tetrahedron_4; akantu::UInt paratype = TETRA1; #endif //AKANTU_USE_IOHELPER akantu::Mesh mesh(3); akantu::MeshIOMSH mesh_io; mesh_io.read("cube1.msh", mesh); akantu::SolidMechanicsModel * model = new akantu::SolidMechanicsModel(mesh); /// model initialization model->initVectors(); /// initialize the vectors akantu::UInt nb_nodes = model->getFEM().getMesh().getNbNodes(); memset(model->getForce().values, 0, 3*nb_nodes*sizeof(akantu::Real)); memset(model->getVelocity().values, 0, 3*nb_nodes*sizeof(akantu::Real)); memset(model->getAcceleration().values, 0, 3*nb_nodes*sizeof(akantu::Real)); memset(model->getDisplacement().values, 0, 3*nb_nodes*sizeof(akantu::Real)); model->readMaterials("material.dat"); model->initMaterials(); model->initModel(); akantu::Real time_step = model->getStableTimeStep(); model->setTimeStep(time_step/10.); model->assembleMassLumped(); std::cout << *model << std::endl; /// boundary conditions akantu::Real eps = 1e-16; for (akantu::UInt i = 0; i < nb_nodes; ++i) { model->getDisplacement().values[3*i] = model->getFEM().getMesh().getNodes().values[3*i] / 100.; if(model->getFEM().getMesh().getNodes().values[3*i] <= eps) { model->getBoundary().values[3*i ] = true; } if(model->getFEM().getMesh().getNodes().values[3*i + 1] <= eps) { model->getBoundary().values[3*i + 1] = true; } } // model->getDisplacement().values[1] = 0.1; #ifdef AKANTU_USE_IOHELPER DumperParaview dumper; // dumper.SetMode(TEXT); dumper.SetPoints(model->getFEM().getMesh().getNodes().values, 3, nb_nodes, "coordinates"); dumper.SetConnectivity((int *)model->getFEM().getMesh().getConnectivity(type).values, paratype, model->getFEM().getMesh().getNbElement(type), C_MODE); dumper.AddNodeDataField(model->getDisplacement().values, 3, "displacements"); dumper.AddNodeDataField(model->getVelocity().values, 3, "velocity"); dumper.AddNodeDataField(model->getMass().values, 1, "mass"); dumper.AddNodeDataField(model->getResidual().values, 3, "force"); dumper.AddElemDataField(model->getMaterial(0).getStrain(type).values, 9, "strain"); dumper.AddElemDataField(model->getMaterial(0).getStress(type).values, 9, "stress"); dumper.SetPrefix("paraview/"); dumper.Init(); #endif //AKANTU_USE_IOHELPER std::ofstream energy; energy.open("energy.csv"); energy << "id,epot,ekin,tot" << std::endl; for(akantu::UInt s = 0; s < max_steps; ++s) { model->explicitPred(); model->updateResidual(); model->updateAcceleration(); model->explicitCorr(); epot = model->getPotentialEnergy(); ekin = model->getKineticEnergy(); std::cerr << "passing step " << s << "/" << max_steps << std::endl; energy << s << "," << epot << "," << ekin << "," << epot + ekin << std::endl; #ifdef AKANTU_USE_IOHELPER if(s % 10 == 0) dumper.Dump(); #endif //AKANTU_USE_IOHELPER } energy.close(); return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_cube3d_tetra10.cc b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_cube3d_tetra10.cc index 206808a53..bc4d886cb 100644 --- a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_cube3d_tetra10.cc +++ b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_cube3d_tetra10.cc @@ -1,139 +1,139 @@ /** * @file test_solid_mechanics_model_cube3d.cc - * @author Guillaume ANCIAUX + * @author Guillaume ANCIAUX * @date Tue Aug 17 11:31:22 2010 * * @brief test of the class SolidMechanicsModel on the 3d cube * * @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 "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "solid_mechanics_model.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER int main(int argc, char *argv[]) { akantu::UInt max_steps = 1000; akantu::Real epot, ekin; #ifdef AKANTU_USE_IOHELPER akantu::ElementType type = akantu::_tetrahedron_10; akantu::UInt paratype = TETRA2; #endif //AKANTU_USE_IOHELPER akantu::Mesh mesh(3); akantu::MeshIOMSH mesh_io; mesh_io.read("cube2.msh", mesh); akantu::SolidMechanicsModel * model = new akantu::SolidMechanicsModel(mesh); /// model initialization model->initVectors(); /// initialize the vectors akantu::UInt nb_nodes = model->getFEM().getMesh().getNbNodes(); memset(model->getForce().values, 0, 3*nb_nodes*sizeof(akantu::Real)); memset(model->getVelocity().values, 0, 3*nb_nodes*sizeof(akantu::Real)); memset(model->getAcceleration().values, 0, 3*nb_nodes*sizeof(akantu::Real)); memset(model->getDisplacement().values, 0, 3*nb_nodes*sizeof(akantu::Real)); model->readMaterials("material.dat"); model->initMaterials(); model->initModel(); akantu::Real time_step = model->getStableTimeStep(); model->setTimeStep(time_step/10.); model->assembleMassLumped(); std::cout << *model << std::endl; /// boundary conditions akantu::Real eps = 1e-2; for (akantu::UInt i = 0; i < nb_nodes; ++i) { model->getDisplacement().values[3*i] = model->getFEM().getMesh().getNodes().values[3*i] / 100.; if(model->getFEM().getMesh().getNodes().values[3*i] <= eps) { model->getBoundary().values[3*i ] = true; } if(model->getFEM().getMesh().getNodes().values[3*i + 1] <= eps) { model->getBoundary().values[3*i + 1] = true; } } // model->getDisplacement().values[1] = 0.1; #ifdef AKANTU_USE_IOHELPER DumperParaview dumper; // dumper.SetMode(TEXT); dumper.SetPoints(model->getFEM().getMesh().getNodes().values, 3, nb_nodes, "coordinates2"); dumper.SetConnectivity((int *)model->getFEM().getMesh().getConnectivity(type).values, paratype, model->getFEM().getMesh().getNbElement(type), C_MODE); dumper.AddNodeDataField(model->getDisplacement().values, 3, "displacements"); dumper.AddNodeDataField(model->getVelocity().values, 3, "velocity"); dumper.AddNodeDataField(model->getMass().values, 1, "mass"); dumper.AddNodeDataField(model->getResidual().values, 3, "force"); dumper.AddElemDataField(model->getMaterial(0).getStrain(type).values, 9, "strain"); dumper.AddElemDataField(model->getMaterial(0).getStress(type).values, 9, "stress"); dumper.SetPrefix("paraview/"); dumper.Init(); #endif //AKANTU_USE_IOHELPER std::ofstream energy; energy.open("energy.csv"); energy << "id,epot,ekin,tot" << std::endl; for(akantu::UInt s = 0; s < max_steps; ++s) { model->explicitPred(); model->updateResidual(); model->updateAcceleration(); model->explicitCorr(); epot = model->getPotentialEnergy(); ekin = model->getKineticEnergy(); std::cerr << "passing step " << s << "/" << max_steps << std::endl; energy << s << "," << epot << "," << ekin << "," << epot + ekin << std::endl; #ifdef AKANTU_USE_IOHELPER if(s % 10 == 0) dumper.Dump(); #endif //AKANTU_USE_IOHELPER } energy.close(); return EXIT_SUCCESS; }