diff --git a/CMakeLists.txt b/CMakeLists.txt index 3d93546ea..4ced2611c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,200 +1,200 @@ #=============================================================================== # @file CMakeLists.txt # @author Nicolas Richart # @date Fri Jun 11 09:46:59 2010 # # @section LICENSE # # # # @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(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 #=============================================================================== # Debug option(AKANTU_DEBUG "Compiles akantu with the debug messages" ON) if(NOT AKANTU_DEBUG) add_definitions(-DAKANTU_NDEBUG) endif(NOT AKANTU_DEBUG) # IOHelper option(AKANTU_USE_IOHELPER "Add IOHelper support in akantu" OFF) if(AKANTU_USE_IOHELPER) include(${AKANTU_CMAKE_DIR}/FindIOHelper.cmake) 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) option(AKANTU_USE_MPI "Add MPI support in akantu" OFF) if(AKANTU_USE_MPI) find_package(MPI) 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) option(AKANTU_USE_SCOTCH "Add Scotch support in akantu" OFF) if(AKANTU_USE_SCOTCH) include(${AKANTU_CMAKE_DIR}/FindScotch.cmake) 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) #=============================================================================== # 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 - common/aka_error.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/material.cc model/solid_mechanics/materials/material_elastic.cc mesh_utils/mesh_io.cc mesh_utils/mesh_io/mesh_io_msh.cc mesh_utils/mesh_partition.cc mesh_utils/mesh_utils.cc synchronizer/ghost_synchronizer.cc synchronizer/synchronizer.cc synchronizer/communicator.cc synchronizer/static_communicator.cc ) -#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_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) set(AKANTU_INCLUDE_DIRS common fem/ mesh_utils/ mesh_utils/mesh_io/ mesh_utils/mesh_partition/ model/ model/integeration_scheme model/solid_mechanics model/solid_mechanics/materials synchronizer/ ) include_directories( ${AKANTU_INCLUDE_DIRS} ${AKANTU_EXTERNAL_LIB_INCLUDE_PATH} ) add_library(akantu ${AKANTU_COMMON_SRC}) set_target_properties(akantu PROPERTIES ${AKANTU_LIBRARY_PROPERTIES}) #=============================================================================== # Tests #=============================================================================== option(AKANTU_TESTS "Activate tests" OFF) if(AKANTU_TESTS) subdirs(test) endif(AKANTU_TESTS) diff --git a/cmake/FindScotch.cmake b/cmake/FindScotch.cmake index 2ac6917d5..62e8488cb 100644 --- a/cmake/FindScotch.cmake +++ b/cmake/FindScotch.cmake @@ -1,46 +1,47 @@ #=============================================================================== # @file FindScotch.cmake # @author Nicolas Richart # @date Tue Aug 25 16:53:57 2010 # # @brief The find_package file for Scotch # # @section LICENSE # # # #=============================================================================== #=============================================================================== set(SCOTCH_LIBRARY "NOTFOUND" CACHE INTERNAL "Cleared" FORCE) find_library(SCOTCH_LIBRARY scotch PATHS ${SCOTCH_DIR} PATH_SUFFIXES src/libscotch lib ) find_library(SCOTCH_LIBRARY_ERR scotcherr PATHS ${SCOTCH_DIR} PATH_SUFFIXES src/libscotch lib ) find_path(SCOTCH_INCLUDE_PATH scotch.h PATHS ${SCOTCH_DIR} PATH_SUFFIXES include scotch src/libscotch ) #=============================================================================== mark_as_advanced(SCOTCH_LIBRARY) +mark_as_advanced(SCOTCH_LIBRARY_ERR) mark_as_advanced(SCOTCH_INCLUDE_PATH) set(SCOTCH_LIBRARIES_ALL ${SCOTCH_LIBRARY} ${SCOTCH_LIBRARY_ERR}) set(SCOTCH_LIBRARIES ${SCOTCH_LIBRARIES_ALL} CACHE INTERNAL "Libraries for scotch" FORCE) #=============================================================================== include(FindPackageHandleStandardArgs) find_package_handle_standard_args(SCOTCH DEFAULT_MSG SCOTCH_LIBRARY SCOTCH_LIBRARY_ERR SCOTCH_INCLUDE_PATH) #=============================================================================== if(NOT SCOTCH_FOUND) set(SCOTCH_DIR "" CACHE PATH "Location of IOHelper source directory.") endif(NOT SCOTCH_FOUND) diff --git a/common/aka_common.cc b/common/aka_common.cc index c05efb523..f87a1bef2 100644 --- a/common/aka_common.cc +++ b/common/aka_common.cc @@ -1,48 +1,49 @@ /** * @file aka_common.cc * @author Nicolas Richart * @date Fri Jun 11 16:56:43 2010 * * @brief Initialization of global variables * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_static_memory.hh" #include "static_communicator.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ void initialize(int * argc, char *** argv) { AKANTU_DEBUG_IN(); StaticMemory::getStaticMemory(); StaticCommunicator * comm = StaticCommunicator::getStaticCommunicator(argc, argv); debug::setParallelContext(comm->whoAmI(), comm->getNbProc()); + debug::initSignalHandler(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void finalize() { AKANTU_DEBUG_IN(); if(StaticMemory::isInstantiated()) delete StaticMemory::getStaticMemory(); if(StaticCommunicator::isInstantiated()) { StaticCommunicator *comm = StaticCommunicator::getStaticCommunicator(); comm->barrier(); delete comm; } AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/common/aka_common.hh b/common/aka_common.hh index b1abc1837..472e91f39 100644 --- a/common/aka_common.hh +++ b/common/aka_common.hh @@ -1,217 +1,228 @@ /** * @file aka_common.hh * @author Nicolas Richart * @date Fri Jun 11 09:48:06 2010 * * @section LICENSE * * * * @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; enum ElementType { _not_defined = 0, _line_1 = 1, // implemented _line_2 = 2, _triangle_1 = 3, _triangle_2 = 4, _tetrahedra_1 = 5, _tetrahedra_2 = 6, - _max_element_type + _max_element_type, + _point }; enum MaterialType { _elastic = 0, _max_material_type }; /* -------------------------------------------------------------------------- */ /* Ghosts handling */ /* -------------------------------------------------------------------------- */ typedef ID SynchronizerID; enum CommunicatorType { _communicator_mpi, _communicator_dummy }; enum GhostSynchronizationTag { /// SolidMechanicsModel tags _gst_smm_mass, _gst_smm_residual, + _gst_smm_boundary, + /// Test _gst_test }; enum GhostType { _not_ghost, _ghost, _casper // not used but a real cute ghost }; +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 _line_1 : stream << "line_1" ; break; case _line_2 : stream << "line_2" ; break; case _triangle_1 : stream << "triangle_1" ; break; case _triangle_2 : stream << "triangle_2" ; break; case _tetrahedra_1 : stream << "tetrahedra_1"; break; case _tetrahedra_2 : stream << "tetrahedra_2"; 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 (); /* -------------------------------------------------------------------------- */ /* 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 = ""; } /* -------------------------------------------------------------------------- */ /* * 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) __END_AKANTU__ #endif /* __AKANTU_COMMON_HH__ */ diff --git a/common/aka_error.cc b/common/aka_error.cc index 3d1a48895..64cb6fb32 100644 --- a/common/aka_error.cc +++ b/common/aka_error.cc @@ -1,32 +1,107 @@ /** * @file aka_error.cc * @author Nicolas Richart * @date Sun Sep 5 21:03:37 2010 * - * @brief + * @brief * * @section LICENSE * * * */ +/* -------------------------------------------------------------------------- */ +#include +#include +#include +#include + /* -------------------------------------------------------------------------- */ #include "aka_error.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ namespace debug { + /* ------------------------------------------------------------------------ */ void setDebugLevel(const DebugLevel & level) { _debug_level = level; } + /* ------------------------------------------------------------------------ */ void setParallelContext(int rank, int size) { std::stringstream sstr; sstr << "[" << std::setfill(' ') << std::right << std::setw(3) << (rank + 1) << "/" << size << "] "; _parallel_context = sstr.str(); } + + /* ------------------------------------------------------------------------ */ + void initSignalHandler() { + struct sigaction action; + + action.sa_handler = &printBacktrace; + sigemptyset(&(action.sa_mask)); + action.sa_flags = SA_RESETHAND; + + sigaction(SIGSEGV, &action, NULL); + } + + /* ------------------------------------------------------------------------ */ + static std::string demangle(const char* symbol) { + std::string trace(symbol); + + std::string::size_type begin = trace.find_first_of('(') + 1; + std::string::size_type end = trace.find_last_of('+'); + + if (begin != std::string::npos && end != std::string::npos) { + std::string sub_trace = trace.substr(begin, end - begin); + + int status; + size_t size; + std::string result; + char * demangled_name; + + if ((demangled_name = abi::__cxa_demangle(sub_trace.c_str(), NULL, &size, &status)) != NULL) { + result = demangled_name; + free(demangled_name); + } + + std::stringstream result_sstr; + result_sstr << trace.substr(0, begin) << result << trace.substr(end); + + return result_sstr.str(); + } + + return trace; + } + + + /* ------------------------------------------------------------------------ */ + void printBacktrace(int sig) { + AKANTU_DEBUG_INFO("Caught signal " << sig << "!"); + + void *array[10]; + size_t size; + char **strings; + size_t i; + + size = backtrace (array, 10); + strings = backtrace_symbols (array, size); + + std::cerr << "BACKTRACE : " << size - 1 << " stack frames." < * @date Mon Jun 14 11:43:22 2010 * * @brief error management and internal exceptions * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_ERROR_HH__ #define __AKANTU_ERROR_HH__ /* -------------------------------------------------------------------------- */ #include #include #ifdef AKANTU_USE_MPI #include #endif /* -------------------------------------------------------------------------- */ #include "aka_common.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ enum DebugLevel { dbl0 = 0, dblError = 0, dblAssert = 0, dbl1 = 1, dblException = 1, dblCritical = 1, dbl2 = 2, dblMajor = 2, dbl3 = 3, dblCall = 3, dblSecondary = 3, dblHead = 3, dbl4 = 4, dblWarning = 4, dbl5 = 5, dblInfo = 5, dbl6 = 6, dblIn = 6, dblOut = 6, dbl7 = 7, dbl8 = 8, dblTrace = 8, dbl9 = 9, dblAccessory = 9, dbl10 = 10, dbl100 = 100, dblDump = 100 }; /* -------------------------------------------------------------------------- */ namespace debug { extern std::ostream & _akantu_cout; extern std::ostream & _akantu_debug_cout; extern DebugLevel _debug_level; extern std::string _parallel_context; - + void setDebugLevel(const DebugLevel & level); void setParallelContext(int rank, int size); + + void initSignalHandler(); + void printBacktrace(int sig); }; /* -------------------------------------------------------------------------- */ /// exception class that can be thrown by akantu class Exception : public std::exception { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: //! full constructor Exception(std::string info, std::string file, unsigned int line) : _info(info), _file(file), _line(line) { } //! destructor virtual ~Exception() throw() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual const char* what() const throw() { std::stringstream stream; stream << "akantu::Exception" << " : " << _info << " [" << _file << ":" << _line << "]"; return stream.str().c_str(); } virtual const char* info() const throw() { return _info.c_str(); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// exception description and additionals std::string _info; /// file it is thrown from std::string _file; /// ligne it is thrown from unsigned int _line; }; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #define AKANTU_LOCATION "(" <<__FILE__ << ":" << __func__ << "():" << __LINE__ << ")" #ifndef AKANTU_USE_MPI #define AKANTU_EXIT(status) \ do { \ - int * p = NULL; \ - *p = 2; \ + akantu::debug::printBacktrace(15); \ exit(status); \ } while(0) #else #define AKANTU_EXIT(status) \ do { \ - int * p = NULL; \ - *p = 2; \ + akantu::debug::printBacktrace(15); \ MPI_Abort(MPI_COMM_WORLD, MPI_ERR_UNKNOWN); \ } while(0) #endif /* -------------------------------------------------------------------------- */ #define AKANTU_EXCEPTION(info) \ do { \ AKANTU_DEBUG(akantu::dblError, "!!! " << info); \ std::stringstream s_info; \ s_info << info ; \ Exception ex(s_info.str(), __FILE__, __LINE__ ); \ throw ex; \ } while(0) /* -------------------------------------------------------------------------- */ #ifdef AKANTU_NDEBUG #define AKANTU_DEBUG_TEST(level) (0) #define AKANTU_DEBUG(level,info) #define AKANTU_DEBUG_IN() #define AKANTU_DEBUG_OUT() #define AKANTU_DEBUG_INFO(info) #define AKANTU_DEBUG_WARNING(info) #define AKANTU_DEBUG_TRACE(info) #define AKANTU_DEBUG_ASSERT(test,info) #define AKANTU_DEBUG_ERROR(info) AKANTU_EXCEPTION(info) /* -------------------------------------------------------------------------- */ #else #define AKANTU_DEBUG(level,info) \ ((::akantu::debug::_debug_level >= level) && \ (::akantu::debug::_akantu_debug_cout \ << ::akantu::debug::_parallel_context \ << info << " " \ << AKANTU_LOCATION \ << std::endl)) #define AKANTU_DEBUG_TEST(level) \ (::akantu::debug::_debug_level >= (level)) #define AKANTU_DEBUG_IN() \ AKANTU_DEBUG(::akantu::dblIn , "==> " << __func__ << "()") #define AKANTU_DEBUG_OUT() \ AKANTU_DEBUG(::akantu::dblOut , "<== " << __func__ << "()") #define AKANTU_DEBUG_INFO(info) \ AKANTU_DEBUG(::akantu::dblInfo , "--- " << info) #define AKANTU_DEBUG_WARNING(info) \ AKANTU_DEBUG(::akantu::dblWarning, "??? " << info) #define AKANTU_DEBUG_TRACE(info) \ AKANTU_DEBUG(::akantu::dblTrace , ">>> " << info) #define AKANTU_DEBUG_ASSERT(test,info) \ do { \ if (!(test)) { \ AKANTU_DEBUG(::akantu::dblAssert, "assert [" << #test << "] " \ << "!!! " << info); \ AKANTU_EXIT(EXIT_FAILURE); \ } \ } while(0) #define AKANTU_DEBUG_ERROR(info) \ do { \ AKANTU_DEBUG(::akantu::dblError, "!!! " << info); \ AKANTU_EXIT(EXIT_FAILURE); \ } while(0) #endif // AKANTU_NDEBUG /* -------------------------------------------------------------------------- */ __END_AKANTU__ #endif /* __AKANTU_ERROR_HH__ */ // LocalWords: acessory diff --git a/fem/element_class.cc b/fem/element_class.cc index 482456b35..b6e275b37 100644 --- a/fem/element_class.cc +++ b/fem/element_class.cc @@ -1,108 +1,116 @@ /** * @file element_class.cc * @author Nicolas Richart * @date Tue Jul 20 10:12:44 2010 * * @brief Common part of element_classes * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #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<> 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<_line_1>::nb_nodes_per_element = 2; template<> ElementType ElementClass<_line_1>::p1_element_type = _line_1; template<> UInt ElementClass<_line_1>::nb_quadrature_points = 1; template<> UInt ElementClass<_line_1>::spatial_dimension = 1; template<> UInt ElementClass<_line_1>::nb_facets = 2; -template<> ElementType ElementClass<_line_1>::facet_type = _not_defined; +template<> ElementType ElementClass<_line_1>::facet_type = _point; template<> UInt ElementClass<_line_1>::vec_facet_connectivity[]= {0, 1}; template<> UInt * ElementClass<_line_1>::facet_connectivity[] = {&vec_facet_connectivity[0], &vec_facet_connectivity[1]}; /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_line_2>::nb_nodes_per_element = 3; template<> ElementType ElementClass<_line_2>::p1_element_type = _line_1; template<> UInt ElementClass<_line_2>::nb_quadrature_points = 2; template<> UInt ElementClass<_line_2>::spatial_dimension = 1; template<> UInt ElementClass<_line_2>::nb_facets = 2; -template<> ElementType ElementClass<_line_2>::facet_type = _not_defined; +template<> ElementType ElementClass<_line_2>::facet_type = _point; template<> UInt ElementClass<_line_2>::vec_facet_connectivity[]= {0, 1}; template<> UInt * ElementClass<_line_2>::facet_connectivity[] = {&vec_facet_connectivity[0], &vec_facet_connectivity[1]}; /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_triangle_1>::nb_nodes_per_element = 3; template<> ElementType ElementClass<_triangle_1>::p1_element_type = _triangle_1; template<> UInt ElementClass<_triangle_1>::nb_quadrature_points = 1; template<> UInt ElementClass<_triangle_1>::spatial_dimension = 2; template<> UInt ElementClass<_triangle_1>::nb_facets = 3; template<> ElementType ElementClass<_triangle_1>::facet_type = _line_1; template<> UInt ElementClass<_triangle_1>::vec_facet_connectivity[]= {0, 1, 1, 2, 2, 0}; template<> UInt * ElementClass<_triangle_1>::facet_connectivity[] = {&vec_facet_connectivity[0], &vec_facet_connectivity[2], &vec_facet_connectivity[4]}; /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_triangle_2>::nb_nodes_per_element = 6; template<> ElementType ElementClass<_triangle_2>::p1_element_type = _triangle_1; template<> UInt ElementClass<_triangle_2>::nb_quadrature_points = 2; template<> UInt ElementClass<_triangle_2>::spatial_dimension = 2; template<> UInt ElementClass<_triangle_2>::nb_facets = 3; template<> ElementType ElementClass<_triangle_2>::facet_type = _line_2; template<> UInt ElementClass<_triangle_2>::vec_facet_connectivity[]= {0, 1, 3, 1, 2, 4, 2, 0, 5}; template<> UInt * ElementClass<_triangle_2>::facet_connectivity[] = {&vec_facet_connectivity[0], &vec_facet_connectivity[3], &vec_facet_connectivity[6]}; /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_tetrahedra_1>::nb_nodes_per_element = 4; template<> ElementType ElementClass<_tetrahedra_1>::p1_element_type = _tetrahedra_1; template<> UInt ElementClass<_tetrahedra_1>::nb_quadrature_points = 1; template<> UInt ElementClass<_tetrahedra_1>::spatial_dimension = 3; template<> UInt ElementClass<_tetrahedra_1>::nb_facets = 4; template<> ElementType ElementClass<_tetrahedra_1>::facet_type = _triangle_1; template<> UInt ElementClass<_tetrahedra_1>::vec_facet_connectivity[]= {0, 2, 1, 1, 2, 3, 2, 0, 3, 0, 1, 3}; template<> UInt * ElementClass<_tetrahedra_1>::facet_connectivity[] = {&vec_facet_connectivity[0], &vec_facet_connectivity[3], &vec_facet_connectivity[6], &vec_facet_connectivity[9]}; /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_tetrahedra_2>::nb_nodes_per_element = 10; template<> ElementType ElementClass<_tetrahedra_2>::p1_element_type = _tetrahedra_1; template<> UInt ElementClass<_tetrahedra_2>::nb_quadrature_points = 4; template<> UInt ElementClass<_tetrahedra_2>::spatial_dimension = 3; template<> UInt ElementClass<_tetrahedra_2>::nb_facets = 4; template<> ElementType ElementClass<_tetrahedra_2>::facet_type = _triangle_2; template<> UInt ElementClass<_tetrahedra_2>::vec_facet_connectivity[]= {0, 2, 1, 6, 5, 4, 1, 2, 3, 5, 8, 9, 2, 0, 3, 8, 7, 6, 0, 1, 3, 4, 9, 7}; template<> UInt * ElementClass<_tetrahedra_2>::facet_connectivity[] = {&vec_facet_connectivity[0], &vec_facet_connectivity[3], &vec_facet_connectivity[6], &vec_facet_connectivity[9]}; /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/fem/fem.cc b/fem/fem.cc index 6a1b7ebed..da561c020 100644 --- a/fem/fem.cc +++ b/fem/fem.cc @@ -1,680 +1,682 @@ /** * @file fem.cc * @author Nicolas Richart * @date Fri Jul 16 11:03:02 2010 * * @brief Implementation of the FEM class * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #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->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 size_of_jacobians = FEM::getJacobianSize(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, size_of_shapes)); std::stringstream sstr_shapesd; sstr_shapesd << id << ":" << ghost << "shapes_derivatives:" << type; Vector * shapes_derivatives_tmp = &(alloc(sstr_shapesd.str(), nb_element, size_of_shapesd)); std::stringstream sstr_jacobians; sstr_jacobians << id << ":" << ghost << "jacobians:" << type; Vector * jacobians_tmp = &(alloc(sstr_jacobians.str(), nb_element, size_of_jacobians)); 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 { \ if (need_rotation) { \ Real local_coord[spatial_dimension * nb_nodes_per_element]; \ Real element_coord[element_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::changeDimension(local_coord, \ spatial_dimension, \ nb_nodes_per_element, \ element_coord); \ \ ElementClass::shapeFunctions(element_coord, \ shapes_val, \ shapesd_val, \ jacobians_val); \ shapes_val += size_of_shapes; \ shapesd_val += size_of_shapesd; \ jacobians_val += size_of_jacobians; \ } \ } \ else { \ Real local_coord[element_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 * element_dimension, \ coord + elem_val[offset + id] * element_dimension, \ element_dimension*sizeof(Real)); \ } \ ElementClass::shapeFunctions(local_coord, \ shapes_val, \ shapesd_val, \ jacobians_val); \ shapes_val += size_of_shapes; \ shapesd_val += size_of_shapesd; \ jacobians_val += size_of_jacobians; \ } \ } \ } while(0) /* -------------------------------------------------------------------------- */ bool need_rotation = mesh->getSpatialDimension() != element_dimension; switch(type) { case _line_1 : { COMPUTE_SHAPES(_line_1 ); break; } case _line_2 : { COMPUTE_SHAPES(_line_2 ); break; } case _triangle_1 : { COMPUTE_SHAPES(_triangle_1 ); break; } case _triangle_2 : { COMPUTE_SHAPES(_triangle_2 ); break; } case _tetrahedra_1 : { COMPUTE_SHAPES(_tetrahedra_1); break; } case _tetrahedra_2 : { COMPUTE_SHAPES(_tetrahedra_2); break; } + case _point: case _not_defined: case _max_element_type: { AKANTU_DEBUG_ERROR("Wrong type : " << type); break; } } #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(); 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 = ""; 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_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 _line_1 : { COMPUTE_NORMALS_ON_QUAD(_line_1 ); break; } case _line_2 : { COMPUTE_NORMALS_ON_QUAD(_line_2 ); break; } case _triangle_1 : { COMPUTE_NORMALS_ON_QUAD(_triangle_1 ); break; } case _triangle_2 : { COMPUTE_NORMALS_ON_QUAD(_triangle_2 ); break; } case _tetrahedra_1 : { COMPUTE_NORMALS_ON_QUAD(_tetrahedra_1); break; } case _tetrahedra_2 : { COMPUTE_NORMALS_ON_QUAD(_tetrahedra_2); break; } + case _point: case _not_defined: case _max_element_type: { AKANTU_DEBUG_ERROR("Wrong type : " << type); break; } } #undef COMPUTE_SHAPES 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 * nb_quadrature_points, "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); Real * shape_val = shapes_loc->values; Real * u_val = in_u.values; Real * uq_val = out_uq.values; UInt offset_uq = out_uq.getNbComponent(); 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; 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; } } 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 * nb_quadrature_points * element_dimension, "The vector out_nablauq(" << out_nablauq.getID() << ") has not the good number of component."); out_nablauq.resize(nb_element); 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; 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 = jacobians[type]; + jac_loc = ghost_jacobians[type]; nb_element = mesh->getNbGhostElement(type); } UInt nb_quadrature_points = FEM::getNbQuadraturePoints(type); UInt size_of_jacobians = FEM::getJacobianSize(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() == mesh->getNbElement(type), - "The vector in_f(" << in_f.getID() - << ") has not the good size."); + AKANTU_DEBUG_ASSERT(in_f.getSize() == nb_element, + "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 * nb_quadrature_points, "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(); 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] * size_of_jacobians; } 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 += size_of_jacobians; } } 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 = shapes_derivatives[type]; + jac_loc = jacobians[type]; nb_element = mesh->getNbElement(type); } else { - jac_loc = ghost_shapes_derivatives[type]; + jac_loc = ghost_jacobians[type]; nb_element = mesh->getNbGhostElement(type); } UInt nb_quadrature_points = FEM::getNbQuadraturePoints(type); UInt size_of_jacobians = FEM::getJacobianSize(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() == mesh->getNbElement(type), + AKANTU_DEBUG_ASSERT(in_f.getSize() == nb_element, "The vector in_f(" << in_f.getID() << ") has not the good size."); AKANTU_DEBUG_ASSERT(in_f.getNbComponent() == nb_quadrature_points, "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(); Real * jac = jac_val; for (UInt el = 0; el < nb_element; ++el) { if(filter_elements != NULL) { jac = jac_val + filter_elem_val[el] * size_of_jacobians; } Real el_intf; 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 += size_of_jacobians; } } 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; 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_inline_impl.cc b/fem/fem_inline_impl.cc index 97bffb622..c1469f9b3 100644 --- a/fem/fem_inline_impl.cc +++ b/fem/fem_inline_impl.cc @@ -1,201 +1,206 @@ /** * @file fem_inline.impl.cc * @author Nicolas Richart * @date Mon Jul 19 12:21:36 2010 * * @brief Implementation of the inline functions of the FEM Class * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ inline void FEM::integrate(const Vector & f, Real * intf, UInt nb_degre_of_freedom, const ElementType & type, const UInt elem, GhostType ghost_type) const { Vector * jac_loc; if(ghost_type == _not_ghost) { jac_loc = shapes_derivatives[type]; } else { jac_loc = ghost_shapes_derivatives[type]; } UInt nb_quadrature_points = FEM::getNbQuadraturePoints(type); UInt size_of_jacobians = FEM::getJacobianSize(type); AKANTU_DEBUG_ASSERT(f.getNbComponent() == nb_degre_of_freedom * nb_quadrature_points, "The vector f do not have the good number of component."); Real * f_val = f.values + elem * f.getNbComponent(); Real * jac_val = jac_loc->values + elem * size_of_jacobians; 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; #define GET_NB_QUAD_POINTS(type) \ nb_quadrature_points = ElementClass::getNbQuadraturePoints() switch(type) { case _line_1 : { GET_NB_QUAD_POINTS(_line_1 ); break; } case _line_2 : { GET_NB_QUAD_POINTS(_line_2 ); break; } case _triangle_1 : { GET_NB_QUAD_POINTS(_triangle_1 ); break; } case _triangle_2 : { GET_NB_QUAD_POINTS(_triangle_2 ); break; } case _tetrahedra_1 : { GET_NB_QUAD_POINTS(_tetrahedra_1); break; } case _tetrahedra_2 : { GET_NB_QUAD_POINTS(_tetrahedra_2); break; } + case _point: case _not_defined: case _max_element_type: { AKANTU_DEBUG_ERROR("Wrong type : " << type); break; } } #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; #define GET_SHAPE_SIZE(type) \ shape_size = ElementClass::getShapeSize() switch(type) { case _line_1 : { GET_SHAPE_SIZE(_line_1 ); break; } case _line_2 : { GET_SHAPE_SIZE(_line_2 ); break; } case _triangle_1 : { GET_SHAPE_SIZE(_triangle_1 ); break; } case _triangle_2 : { GET_SHAPE_SIZE(_triangle_2 ); break; } case _tetrahedra_1 : { GET_SHAPE_SIZE(_tetrahedra_1); break; } case _tetrahedra_2 : { GET_SHAPE_SIZE(_tetrahedra_2); break; } + case _point: case _not_defined: case _max_element_type: { AKANTU_DEBUG_ERROR("Wrong type : " << type); break; } } #undef GET_SHAPE_SIZE AKANTU_DEBUG_OUT(); return shape_size; } /* -------------------------------------------------------------------------- */ inline UInt FEM::getShapeDerivativesSize(const ElementType & type) { AKANTU_DEBUG_IN(); UInt shape_derivatives_size; #define GET_SHAPE_DERIVATIVES_SIZE(type) \ shape_derivatives_size = ElementClass::getShapeDerivativesSize() switch(type) { case _line_1 : { GET_SHAPE_DERIVATIVES_SIZE(_line_1 ); break; } case _line_2 : { GET_SHAPE_DERIVATIVES_SIZE(_line_2 ); break; } case _triangle_1 : { GET_SHAPE_DERIVATIVES_SIZE(_triangle_1 ); break; } case _triangle_2 : { GET_SHAPE_DERIVATIVES_SIZE(_triangle_2 ); break; } case _tetrahedra_1 : { GET_SHAPE_DERIVATIVES_SIZE(_tetrahedra_1); break; } case _tetrahedra_2 : { GET_SHAPE_DERIVATIVES_SIZE(_tetrahedra_2); break; } + case _point: case _not_defined: case _max_element_type: { AKANTU_DEBUG_ERROR("Wrong type : " << type); break; } } #undef GET_SHAPE_DERIVATIVES_SIZE AKANTU_DEBUG_OUT(); return shape_derivatives_size; } /* -------------------------------------------------------------------------- */ inline UInt FEM::getJacobianSize(const ElementType & type) { AKANTU_DEBUG_IN(); UInt jacobian_size; #define GET_JACOBIAN_SIZE(type) \ jacobian_size = ElementClass::getJacobianSize() switch(type) { case _line_1 : { GET_JACOBIAN_SIZE(_line_1 ); break; } case _line_2 : { GET_JACOBIAN_SIZE(_line_2 ); break; } case _triangle_1 : { GET_JACOBIAN_SIZE(_triangle_1 ); break; } case _triangle_2 : { GET_JACOBIAN_SIZE(_triangle_2 ); break; } case _tetrahedra_1 : { GET_JACOBIAN_SIZE(_tetrahedra_1); break; } case _tetrahedra_2 : { GET_JACOBIAN_SIZE(_tetrahedra_2); break; } + case _point: case _not_defined: case _max_element_type: { AKANTU_DEBUG_ERROR("Wrong type : " << type); break; } } #undef GET_JACOBIAN_SIZE AKANTU_DEBUG_OUT(); return jacobian_size; } /* -------------------------------------------------------------------------- */ inline Real FEM::getElementInradius(Real * coord, const ElementType & type) { AKANTU_DEBUG_IN(); Real inradius; #define GET_INRADIUS(type) \ inradius = ElementClass::getInradius(coord); \ switch(type) { case _line_1 : { GET_INRADIUS(_line_1 ); break; } case _line_2 : { GET_INRADIUS(_line_2 ); break; } case _triangle_1 : { GET_INRADIUS(_triangle_1 ); break; } case _triangle_2 : { GET_INRADIUS(_triangle_2 ); break; } case _tetrahedra_1 : { GET_INRADIUS(_tetrahedra_1); break; } case _tetrahedra_2 : { GET_INRADIUS(_tetrahedra_2); break; } + case _point: case _not_defined: case _max_element_type: { AKANTU_DEBUG_ERROR("Wrong type : " << type); break; } } #undef GET_INRADIUS AKANTU_DEBUG_OUT(); return inradius; } /* -------------------------------------------------------------------------- */ diff --git a/fem/mesh_inline_impl.cc b/fem/mesh_inline_impl.cc index 978e44761..b116a2fea 100644 --- a/fem/mesh_inline_impl.cc +++ b/fem/mesh_inline_impl.cc @@ -1,385 +1,391 @@ /** * @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 * * * */ /* -------------------------------------------------------------------------- */ 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::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 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 UInt Mesh::getNbNodesPerElement(const ElementType & type) { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element; #define GET_NB_NODES_PER_ELEMENT(type) \ nb_nodes_per_element = ElementClass::getNbNodesPerElement() switch(type) { case _line_1 : { GET_NB_NODES_PER_ELEMENT(_line_1 ); break; } case _line_2 : { GET_NB_NODES_PER_ELEMENT(_line_2 ); break; } case _triangle_1 : { GET_NB_NODES_PER_ELEMENT(_triangle_1 ); break; } case _triangle_2 : { GET_NB_NODES_PER_ELEMENT(_triangle_2 ); break; } case _tetrahedra_1 : { GET_NB_NODES_PER_ELEMENT(_tetrahedra_1); break; } case _tetrahedra_2 : { GET_NB_NODES_PER_ELEMENT(_tetrahedra_2); break; } + case _point : { GET_NB_NODES_PER_ELEMENT(_point ); break; } case _not_defined: case _max_element_type: { AKANTU_DEBUG_ERROR("Wrong type : " << type); break; } } #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; #define GET_ELEMENT_P1(type) \ element_p1 = ElementClass::getP1ElementType() switch(type) { case _line_1 : { GET_ELEMENT_P1(_line_1 ); break; } case _line_2 : { GET_ELEMENT_P1(_line_2 ); break; } case _triangle_1 : { GET_ELEMENT_P1(_triangle_1 ); break; } case _triangle_2 : { GET_ELEMENT_P1(_triangle_2 ); break; } case _tetrahedra_1 : { GET_ELEMENT_P1(_tetrahedra_1); break; } case _tetrahedra_2 : { GET_ELEMENT_P1(_tetrahedra_2); break; } + case _point : { GET_ELEMENT_P1(_point ); break; } case _not_defined: case _max_element_type: { AKANTU_DEBUG_ERROR("Wrong type : " << type); break; } } #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; #define GET_SPATIAL_DIMENSION(type) \ spatial_dimension = ElementClass::getSpatialDimension() switch(type) { case _line_1 : { GET_SPATIAL_DIMENSION(_line_1 ); break; } case _line_2 : { GET_SPATIAL_DIMENSION(_line_2 ); break; } case _triangle_1 : { GET_SPATIAL_DIMENSION(_triangle_1 ); break; } case _triangle_2 : { GET_SPATIAL_DIMENSION(_triangle_2 ); break; } case _tetrahedra_1 : { GET_SPATIAL_DIMENSION(_tetrahedra_1); break; } case _tetrahedra_2 : { GET_SPATIAL_DIMENSION(_tetrahedra_2); break; } + case _point : { GET_SPATIAL_DIMENSION(_point ); break; } case _not_defined: case _max_element_type: { AKANTU_DEBUG_ERROR("Wrong type : " << type); break; } } #undef GET_SPATIAL_DIMENSION AKANTU_DEBUG_OUT(); return spatial_dimension; } /* -------------------------------------------------------------------------- */ inline const ElementType Mesh::getFacetElementType(const ElementType & type) { AKANTU_DEBUG_IN(); ElementType surface_type; #define GET_FACET_TYPE(type) \ surface_type = ElementClass::getFacetElementType() switch(type) { case _line_1 : { GET_FACET_TYPE(_line_1 ); break; } case _line_2 : { GET_FACET_TYPE(_line_2 ); break; } case _triangle_1 : { GET_FACET_TYPE(_triangle_1 ); break; } case _triangle_2 : { GET_FACET_TYPE(_triangle_2 ); break; } case _tetrahedra_1 : { GET_FACET_TYPE(_tetrahedra_1); break; } case _tetrahedra_2 : { GET_FACET_TYPE(_tetrahedra_2); break; } + case _point : { GET_FACET_TYPE(_point ); break; } case _not_defined: case _max_element_type: { AKANTU_DEBUG_ERROR("Wrong type : " << type); break; } } -#undef GET_SURFACE_TYPE +#undef GET_FACET_TYPE AKANTU_DEBUG_OUT(); return surface_type; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbFacetsPerElementType(const ElementType & type) const { AKANTU_DEBUG_IN(); UInt n_facet; #define GET_NB_FACET(type) \ n_facet = ElementClass::getNbFacetsPerElement() switch(type) { case _line_1 : { GET_NB_FACET(_line_1 ); break; } case _line_2 : { GET_NB_FACET(_line_2 ); break; } case _triangle_1 : { GET_NB_FACET(_triangle_1 ); break; } case _triangle_2 : { GET_NB_FACET(_triangle_2 ); break; } case _tetrahedra_1 : { GET_NB_FACET(_tetrahedra_1); break; } case _tetrahedra_2 : { GET_NB_FACET(_tetrahedra_2); break; } + case _point : { GET_NB_FACET(_point ); break; } case _not_defined: case _max_element_type: { AKANTU_DEBUG_ERROR("Wrong type : " << type); break; } } -#undef GET_SURFACE_TYPE +#undef GET_NB_FACET AKANTU_DEBUG_OUT(); return n_facet; } /* -------------------------------------------------------------------------- */ inline UInt ** Mesh::getFacetLocalConnectivityPerElementType(const ElementType & type) const { AKANTU_DEBUG_IN(); UInt ** facet_conn; #define GET_FACET_CON(type) \ facet_conn = ElementClass::getFacetLocalConnectivityPerElement() switch(type) { case _line_1 : { GET_FACET_CON(_line_1 ); break; } case _line_2 : { GET_FACET_CON(_line_2 ); break; } case _triangle_1 : { GET_FACET_CON(_triangle_1 ); break; } case _triangle_2 : { GET_FACET_CON(_triangle_2 ); break; } case _tetrahedra_1 : { GET_FACET_CON(_tetrahedra_1); break; } case _tetrahedra_2 : { GET_FACET_CON(_tetrahedra_2); break; } + case _point : { GET_FACET_CON(_point ); break; } case _not_defined: case _max_element_type: { AKANTU_DEBUG_ERROR("Wrong type : " << type); break; } } -#undef GET_SURFACE_TYPE +#undef GET_FACET_CON AKANTU_DEBUG_OUT(); return facet_conn; } diff --git a/mesh_utils/mesh_utils.cc b/mesh_utils/mesh_utils.cc index 5129e3d51..0c24a4c44 100644 --- a/mesh_utils/mesh_utils.cc +++ b/mesh_utils/mesh_utils.cc @@ -1,321 +1,321 @@ /** * @file mesh_utils.cc * @author Guillaume ANCIAUX * @date Wed Aug 18 14:21:00 2010 * * @brief All mesh utils necessary for various tasks * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #include "mesh_utils.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ void MeshUtils::buildNode2Elements(const Mesh & mesh, Vector & node_offset, Vector & node_to_elem) { AKANTU_DEBUG_IN(); 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) != mesh.getSpatialDimension()) 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 for (UInt t = 0; t < nb_good_types; ++t) { memset(node_offset_val, 0, (nb_nodes + 1)*sizeof(UInt)); 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::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; + // 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.getNbFacetsPerElementType(type); node_in_facet[nb_good_types] = mesh.getFacetLocalConnectivityPerElementType(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 happend"); //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(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt * conn = local_connectivities; // Vector nodes_numbers; // UInt node_counter = 0; 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;//node_counter; // node_counter++; } else { *conn = nn; } *conn++; } } 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)); // if(old_nodes) { // old_nodes->resize(nodes_numbers.getSize()); // // std::map::iterator nn; // for(nn = nodes_numbers.begin(); // nn != nodes_numbers.end(); // ++nn) { // old_nodes->push_back(nn->first); // } // } AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/model/model.hh b/model/model.hh index e697baa5f..85408ee9c 100644 --- a/model/model.hh +++ b/model/model.hh @@ -1,96 +1,98 @@ /** * @file model.hh * @author Nicolas Richart * @date Thu Jul 22 11:02:42 2010 * * @brief Interface of a model * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MODEL_HH__ #define __AKANTU_MODEL_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "mesh.hh" #include "fem.hh" #include "mesh_utils.hh" #include "ghost_synchronizer.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ class Model : public Memory, public GhostSynchronizer { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: inline Model(Mesh & mesh, UInt spatial_dimension = 0, const ModelID & id = "model", const MemoryID & memory_id = 0); inline virtual ~Model(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual void initModel() = 0; /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const = 0; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); AKANTU_GET_MACRO(FEM, *fem, const FEM &); inline FEM & getFEMBoundary(); + virtual void synchronizeBoundaries() {}; + /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id ModelID id; /// the spatial dimension UInt spatial_dimension; /// the main fem object present in all models FEM * fem; /// the fem object present in all models for boundaries FEM * fem_boundary; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "model_inline_impl.cc" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Model & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_MODEL_HH__ */ diff --git a/model/solid_mechanics/material.cc b/model/solid_mechanics/material.cc index 81c04202b..3f4ca5741 100644 --- a/model/solid_mechanics/material.cc +++ b/model/solid_mechanics/material.cc @@ -1,199 +1,220 @@ /** * @file material.cc * @author Nicolas Richart * @date Tue Jul 27 11:43:41 2010 * * @brief Implementation of the common part of the material class * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "solid_mechanics_model.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ Material::Material(SolidMechanicsModel & model, const MaterialID & id) : Memory(model.getMemoryID()), id(id), spatial_dimension(model.getSpatialDimension()), name(""), model(&model), potential_energy_flag(false), potential_energy_vector(false), is_init(false) { AKANTU_DEBUG_IN(); for(UInt t = _not_defined; t < _max_element_type; ++t) { this->stress [t] = NULL; this->strain [t] = NULL; this->potential_energy [t] = NULL; this->element_filter [t] = NULL; this->ghost_stress [t] = NULL; this->ghost_strain [t] = NULL; this->ghost_potential_energy[t] = NULL; this->ghost_element_filter [t] = NULL; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Material::~Material() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::setParam(const std::string & key, const std::string & value, __attribute__ ((unused)) const MaterialID & id) { if(key == "name") name = std::string(value); else AKANTU_DEBUG_ERROR("Unknown material property : " << key); } /* -------------------------------------------------------------------------- */ void Material::initMaterial() { AKANTU_DEBUG_IN(); /// for each connectivity types allocate the element filer array of the material UInt spatial_dimension = model->getSpatialDimension(); const Mesh::ConnectivityTypeList & type_list = model->getFEM().getMesh().getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; std::stringstream sstr; sstr << id << ":element_filer:"<< *it; element_filter[*it] = &(alloc (sstr.str(), 0, 1)); UInt nb_quadrature_points = FEM::getNbQuadraturePoints(*it); std::stringstream sstr_stre; sstr_stre << id << ":stress:" << *it; std::stringstream sstr_stra; sstr_stra << id << ":strain:" << *it; stress[*it] = &(alloc(sstr_stre.str(), 0, nb_quadrature_points * spatial_dimension * spatial_dimension)); strain[*it] = &(alloc(sstr_stra.str(), 0, nb_quadrature_points * spatial_dimension * spatial_dimension)); } const Mesh::ConnectivityTypeList & ghost_type_list = model->getFEM().getMesh().getConnectivityTypeList(_ghost); for(it = ghost_type_list.begin(); it != ghost_type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; std::stringstream sstr; sstr << id << ":ghost_element_filer:"<< *it; ghost_element_filter[*it] = &(alloc (sstr.str(), 0, 1)); UInt nb_quadrature_points = FEM::getNbQuadraturePoints(*it); std::stringstream sstr_stre; sstr_stre << id << ":ghost_stress:" << *it; std::stringstream sstr_stra; sstr_stra << id << ":ghost_strain:" << *it; ghost_stress[*it] = &(alloc(sstr_stre.str(), 0, nb_quadrature_points * spatial_dimension * spatial_dimension)); ghost_strain[*it] = &(alloc(sstr_stra.str(), 0, nb_quadrature_points * spatial_dimension * spatial_dimension)); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::updateResidual(Vector & current_position, GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Vector & residual = model->getResidual(); const Mesh::ConnectivityTypeList & type_list = model->getFEM().getMesh().getConnectivityTypeList(ghost_type); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); UInt size_of_shapes_derivatives = FEM::getShapeDerivativesSize(*it); UInt nb_quadrature_points = FEM::getNbQuadraturePoints(*it); Vector * strain_vect; Vector * stress_vect; Vector * elem_filter; const Vector * shapes_derivatives; if(ghost_type == _not_ghost) { elem_filter = element_filter[*it]; strain_vect = strain[*it]; stress_vect = stress[*it]; shapes_derivatives = &(model->getFEM().getShapesDerivatives(*it)); } else { elem_filter = ghost_element_filter[*it]; strain_vect = ghost_strain[*it]; stress_vect = ghost_stress[*it]; shapes_derivatives = &(model->getFEM().getGhostShapesDerivatives(*it)); } UInt nb_element = elem_filter->getSize(); /// compute the deformation gradient (current_position -> strain) model->getFEM().gradientOnQuadraturePoints(current_position, *strain_vect, spatial_dimension, *it, ghost_type, elem_filter); /// (strain -> stress) constitutiveLaw(*it, ghost_type); /// compute \sigma * \partial \phi / \partial X Vector * sigma_dphi_dx = new Vector(nb_element, size_of_shapes_derivatives, "sigma_x_dphi_/_dX"); Real * shapesd = shapes_derivatives->values; UInt size_of_shapesd = shapes_derivatives->getNbComponent(); Real * shapesd_val; Real * stress_val = stress_vect->values; Real * sigma_dphi_dx_val = sigma_dphi_dx->values; UInt * elem_filter_val = elem_filter->values; UInt offset_shapesd_val = spatial_dimension * nb_nodes_per_element; UInt offset_stress_val = spatial_dimension * spatial_dimension; UInt offset_sigma_dphi_dx_val = spatial_dimension * nb_nodes_per_element; for (UInt el = 0; el < nb_element; ++el) { shapesd_val = shapesd + elem_filter_val[el]*size_of_shapesd; for (UInt q = 0; q < nb_quadrature_points; ++q) { Math::matrix_matrixt(nb_nodes_per_element, spatial_dimension, spatial_dimension, shapesd_val, stress_val, sigma_dphi_dx_val); shapesd_val += offset_shapesd_val; stress_val += offset_stress_val; sigma_dphi_dx_val += offset_sigma_dphi_dx_val; } } /// integrate \sigma * \partial \phi / \partial X Vector * int_sigma_dphi_dx = new Vector(nb_element, nb_nodes_per_element * spatial_dimension, "int_sigma_x_dphi_/_dX"); model->getFEM().integrate(*sigma_dphi_dx, *int_sigma_dphi_dx, size_of_shapes_derivatives / nb_quadrature_points, - *it, ghost_type, element_filter[*it]); + *it, ghost_type, + elem_filter); delete sigma_dphi_dx; /// assemble model->getFEM().assembleVector(*int_sigma_dphi_dx, residual, residual.getNbComponent(), - *it, ghost_type, element_filter[*it], -1); + *it, ghost_type, elem_filter, -1); delete int_sigma_dphi_dx; } AKANTU_DEBUG_OUT(); } + +/* -------------------------------------------------------------------------- */ +Real Material::getPotentialEnergy() { + AKANTU_DEBUG_IN(); + Real epot = 0.; + + /// fill the element filters of the materials using the element_material arrays + const Mesh::ConnectivityTypeList & type_list = model->getFEM().getMesh().getConnectivityTypeList(_not_ghost); + Mesh::ConnectivityTypeList::const_iterator it; + for(it = type_list.begin(); it != type_list.end(); ++it) { + if(model->getFEM().getMesh().getSpatialDimension(*it) != spatial_dimension) continue; + + epot += model->getFEM().integrate(*potential_energy[*it], *it, + _not_ghost, element_filter[*it]); + } + + AKANTU_DEBUG_OUT(); + return epot; +} + /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/model/solid_mechanics/material.hh b/model/solid_mechanics/material.hh index 0d1fc6260..77cc33b27 100644 --- a/model/solid_mechanics/material.hh +++ b/model/solid_mechanics/material.hh @@ -1,178 +1,238 @@ /** * @file material.hh * @author Nicolas Richart * @date Fri Jul 23 09:06:29 2010 * * @brief Mother class for all materials * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_HH__ #define __AKANTU_MATERIAL_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "fem.hh" #include "mesh.hh" //#include "solid_mechanics_model.hh" - +#include "static_communicator.hh" /* -------------------------------------------------------------------------- */ namespace akantu { class SolidMechanicsModel; }; __BEGIN_AKANTU__ class Material : public Memory { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Material(SolidMechanicsModel & model, const MaterialID & id = ""); virtual ~Material(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// read properties virtual void setParam(const std::string & key, const std::string & value, const MaterialID & id); /// initialize the material computed parameter virtual void initMaterial(); /// compute the residual for this material void updateResidual(Vector & current_position, GhostType ghost_type = _not_ghost); /// constitutive law virtual void constitutiveLaw(ElementType el_type, GhostType ghost_type = _not_ghost) = 0; /// compute the stable time step for an element of size h virtual Real getStableTimeStep(Real h) = 0; /// add an element to the local mesh filter inline void addElement(ElementType type, UInt element); /// add an element to the local mesh filter for ghost element inline void addGhostElement(ElementType type, UInt element); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const = 0; /* ------------------------------------------------------------------------ */ /* Ghost Synchronizer inherited members */ /* ------------------------------------------------------------------------ */ public: inline virtual UInt getNbDataToPack(const Element & element, GhostSynchronizationTag tag); inline virtual UInt getNbDataToUnpack(const Element & element, GhostSynchronizationTag tag); inline virtual void packData(Real ** buffer, const Element & element, GhostSynchronizationTag tag); inline virtual void unpackData(Real ** buffer, const Element & element, GhostSynchronizationTag tag); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(ID, id, const MaterialID &); AKANTU_GET_MACRO(Rho, rho, Real); void setPotentialEnergyFlagOn(); void setPotentialEnergyFlagOff(); + Real getPotentialEnergy(); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Strain, strain, const Vector &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Stress, stress, const Vector &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(PotentialEnergy, potential_energy, const Vector &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id of the material MaterialID id; /// spatial dimension UInt spatial_dimension; /// material name std::string name; /// The model to witch the material belong SolidMechanicsModel * model; /// density : rho Real rho; /// stresses arrays ordered by element types ByElementTypeReal stress; /// strains arrays ordered by element types ByElementTypeReal strain; /// list of element handled by the material ByElementTypeUInt element_filter; /// stresses arrays ordered by element types ByElementTypeReal ghost_stress; /// strains arrays ordered by element types ByElementTypeReal ghost_strain; /// list of element handled by the material ByElementTypeUInt ghost_element_filter; /// has to compute potential energy or not bool potential_energy_flag; /// is the vector for potential energy initialized bool potential_energy_vector; /// potential energy by element ByElementTypeReal potential_energy; /// potential energy by element ByElementTypeReal ghost_potential_energy; /// boolean to know if the material has been initialized bool is_init; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_inline_impl.cc" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Material & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #include "materials/material_elastic.hh" -#endif /* __AKANTU_MATERIAL_HH__ */ +/* -------------------------------------------------------------------------- */ +/* Auto loop */ +/* -------------------------------------------------------------------------- */ +#define MATERIAL_LITTLE_DISP_QUADRATURE_POINT_LOOP_BEGIN \ + UInt nb_quadrature_points = FEM::getNbQuadraturePoints(el_type); \ + UInt size_strain = spatial_dimension * spatial_dimension; \ + \ + UInt nb_element; \ + Real * strain_val; \ + Real * stress_val; \ + bool potential_energy_flag_tmp; \ + \ + if(ghost_type == _not_ghost) { \ + nb_element = element_filter[el_type]->getSize(); \ + strain_val = strain[el_type]->values; \ + stress_val = stress[el_type]->values; \ + } else { \ + nb_element = ghost_element_filter[el_type]->getSize(); \ + strain_val = ghost_strain[el_type]->values; \ + stress_val = ghost_stress[el_type]->values; \ + potential_energy_flag_tmp = potential_energy_flag; \ + potential_energy_flag = false; \ + } \ + \ + if (nb_element == 0) return; \ + \ + Real * epot = NULL; \ + if (potential_energy_flag) epot = potential_energy[el_type]->values; \ + \ + Real F[3*3]; \ + Real sigma[3*3]; \ + \ + for (UInt el = 0; el < nb_element; ++el) { \ + for (UInt q = 0; q < nb_quadrature_points; ++q) { \ + memset(F, 0, 3 * 3 * sizeof(Real)); \ + \ + for (UInt i = 0; i < spatial_dimension; ++i) \ + for (UInt j = 0; j < spatial_dimension; ++j) \ + F[3*i + j] = strain_val[spatial_dimension * i + j]; \ + \ + for (UInt i = 0; i < spatial_dimension; ++i) F[i*3 + i] -= 1; + + +#define MATERIAL_LITTLE_DISP_QUADRATURE_POINT_LOOP_END \ + for (UInt i = 0; i < spatial_dimension; ++i) \ + for (UInt j = 0; j < spatial_dimension; ++j) \ + stress_val[spatial_dimension*i + j] = sigma[3 * i + j]; \ + \ + strain_val += size_strain; \ + stress_val += size_strain; \ + if (potential_energy_flag) epot += nb_quadrature_points; \ + } \ + } \ + \ + if(ghost_type == _ghost) { \ + potential_energy_flag = potential_energy_flag_tmp; \ + } + + +#endif /* __AKANTU_MATERIAL_HH__ */ diff --git a/model/solid_mechanics/material_inline_impl.cc b/model/solid_mechanics/material_inline_impl.cc index 23029729e..09c82bc67 100644 --- a/model/solid_mechanics/material_inline_impl.cc +++ b/model/solid_mechanics/material_inline_impl.cc @@ -1,109 +1,110 @@ /** * @file material_inline_impl.cc * @author Nicolas Richart * @date Tue Jul 27 11:57:43 2010 * * @brief Implementation of the inline functions of the class material * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ inline void Material::addElement(ElementType type, UInt element) { element_filter[type]->push_back(element); strain[type]->push_back(REAL_INIT_VALUE); stress[type]->push_back(REAL_INIT_VALUE); if(potential_energy_vector) potential_energy[type]->push_back(REAL_INIT_VALUE); } /* -------------------------------------------------------------------------- */ inline void Material::addGhostElement(ElementType type, UInt element) { ghost_element_filter[type]->push_back(element); ghost_strain[type]->push_back(REAL_INIT_VALUE); ghost_stress[type]->push_back(REAL_INIT_VALUE); - if(potential_energy_vector) - ghost_potential_energy[type]->push_back(REAL_INIT_VALUE); + + // if(potential_energy_vector) + // ghost_potential_energy[type]->push_back(REAL_INIT_VALUE); } /* -------------------------------------------------------------------------- */ inline void Material::setPotentialEnergyFlagOn(){ AKANTU_DEBUG_IN(); if(!potential_energy_vector) { /// for each connectivity types allocate the element filer array of the material for(UInt t = _not_defined + 1; t < _max_element_type; ++t) { ElementType type = (ElementType) t; if(Mesh::getSpatialDimension(type) != spatial_dimension) continue; UInt nb_quadrature_points = FEM::getNbQuadraturePoints(type); if(element_filter[type] != NULL) { UInt nb_element = element_filter[type]->getSize(); std::stringstream sstr; sstr << id << ":potential_energy:"<< type; potential_energy[type] = &(alloc (sstr.str(), nb_element, nb_quadrature_points, REAL_INIT_VALUE)); } - if(ghost_element_filter[type] != NULL) { - UInt nb_element = ghost_element_filter[type]->getSize(); - std::stringstream sstr; sstr << id << ":ghost_potential_energy:"<< type; - ghost_potential_energy[type] = &(alloc (sstr.str(), nb_element, - nb_quadrature_points, - REAL_INIT_VALUE)); - } + // if(ghost_element_filter[type] != NULL) { + // UInt nb_element = ghost_element_filter[type]->getSize(); + // std::stringstream sstr; sstr << id << ":ghost_potential_energy:"<< type; + // ghost_potential_energy[type] = &(alloc (sstr.str(), nb_element, + // nb_quadrature_points, + // REAL_INIT_VALUE)); + // } } potential_energy_vector = true; } potential_energy_flag = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline void Material::setPotentialEnergyFlagOff(){ AKANTU_DEBUG_IN(); potential_energy_flag = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline UInt Material::getNbDataToPack(const Element & element, GhostSynchronizationTag tag) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); return 0; }; /* -------------------------------------------------------------------------- */ inline UInt Material::getNbDataToUnpack(const Element & element, GhostSynchronizationTag tag) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); return 0; }; /* -------------------------------------------------------------------------- */ inline void Material::packData(Real ** buffer, const Element & element, GhostSynchronizationTag tag) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); }; /* -------------------------------------------------------------------------- */ inline void Material::unpackData(Real ** buffer, const Element & element, GhostSynchronizationTag tag) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); }; diff --git a/model/solid_mechanics/materials/material_elastic.cc b/model/solid_mechanics/materials/material_elastic.cc index abce6720f..c40695505 100644 --- a/model/solid_mechanics/materials/material_elastic.cc +++ b/model/solid_mechanics/materials/material_elastic.cc @@ -1,131 +1,87 @@ /** * @file material_elastic.cc * @author Nicolas Richart * @date Tue Jul 27 11:53:52 2010 * * @brief Specialization of the material class for the elastic material * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #include "material_elastic.hh" #include "solid_mechanics_model.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ MaterialElastic::MaterialElastic(SolidMechanicsModel & model, const MaterialID & id) : Material(model, id) { AKANTU_DEBUG_IN(); rho = 0; E = 0; nu = 1./2.; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialElastic::initMaterial() { AKANTU_DEBUG_IN(); Material::initMaterial(); lambda = nu * E / ((1 + nu) * (1 - 2*nu)); mu = E / (2 * (1 + nu)); kpa = lambda + 2./3. * mu; is_init = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialElastic::constitutiveLaw(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); - UInt nb_quadrature_points = FEM::getNbQuadraturePoints(el_type); - UInt size_strain = spatial_dimension * spatial_dimension; - - UInt nb_element; - Real * strain_val; - Real * stress_val; - - if(ghost_type == _not_ghost) { - nb_element = element_filter[el_type]->getSize(); - strain_val = strain[el_type]->values; - stress_val = stress[el_type]->values; - } else { - nb_element = ghost_element_filter[el_type]->getSize(); - strain_val = ghost_strain[el_type]->values; - stress_val = ghost_stress[el_type]->values; - potential_energy_flag = false; - } - - if (nb_element == 0) return; - - Real * epot = NULL; - if (potential_energy_flag) epot = potential_energy[el_type]->values; - - Real F[3*3]; - Real sigma[3*3]; - - /// for each element - for (UInt el = 0; el < nb_element; ++el) { - /// for each quadrature points - for (UInt q = 0; q < nb_quadrature_points; ++q) { - memset(F, 0, 3 * 3 * sizeof(Real)); - for (UInt i = 0; i < spatial_dimension; ++i) - for (UInt j = 0; j < spatial_dimension; ++j) - F[3*i + j] = strain_val[spatial_dimension * i + j]; + MATERIAL_LITTLE_DISP_QUADRATURE_POINT_LOOP_BEGIN; + constitutiveLaw(F, sigma, epot); + MATERIAL_LITTLE_DISP_QUADRATURE_POINT_LOOP_END; - for (UInt i = 0; i < spatial_dimension; ++i) F[i*3 + i] -= 1; - - constitutiveLaw(F, sigma, epot); - - for (UInt i = 0; i < spatial_dimension; ++i) - for (UInt j = 0; j < spatial_dimension; ++j) - stress_val[spatial_dimension*i + j] = sigma[3 * i + j]; - - strain_val += size_strain; - stress_val += size_strain; - if (potential_energy_flag) epot += nb_quadrature_points; - } - } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialElastic::setParam(const std::string & key, const std::string & value, const MaterialID & id) { std::stringstream sstr(value); if(key == "rho") { sstr >> rho; } else if(key == "E") { sstr >> E; } else if(key == "nu") { sstr >> nu; } else { Material::setParam(key, value, id); } } /* -------------------------------------------------------------------------- */ void MaterialElastic::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "Material<_elastic> [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + name : " << name << std::endl; stream << space << " + density : " << rho << std::endl; stream << space << " + Young's modulus : " << E << std::endl; stream << space << " + Poisson's ratio : " << nu << std::endl; if(is_init) { stream << space << " + First Lamé coefficient : " << lambda << std::endl; stream << space << " + Second Lamé coefficient : " << mu << std::endl; stream << space << " + Bulk coefficient : " << kpa << std::endl; } stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/model/solid_mechanics/solid_mechanics_model.cc b/model/solid_mechanics/solid_mechanics_model.cc index 82aa33964..7e0c0f6fa 100644 --- a/model/solid_mechanics/solid_mechanics_model.cc +++ b/model/solid_mechanics/solid_mechanics_model.cc @@ -1,575 +1,618 @@ /** * @file solid_mechanics_model.cc * @author Nicolas Richart * @date Thu Jul 22 14:35:38 2010 * * @brief Implementation of the SolidMechanicsModel class * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" #include "material.hh" #include "aka_math.hh" #include "integration_scheme/central_difference.hh" +#include "static_communicator.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh, UInt spatial_dimension, const ModelID & id, const MemoryID & memory_id) : Model(mesh, spatial_dimension, id, memory_id), time_step(NAN), f_m2a(1.0), integrator(new CentralDifference()) { AKANTU_DEBUG_IN(); this->displacement = NULL; this->mass = NULL; this->velocity = NULL; this->acceleration = NULL; this->force = NULL; this->residual = NULL; this->boundary = NULL; for(UInt t = _not_defined; t < _max_element_type; ++t) { this->element_material[t] = NULL; this->ghost_element_material[t] = NULL; } registerTag(_gst_smm_mass, "Mass"); registerTag(_gst_smm_residual, "Explicit Residual"); + registerTag(_gst_smm_boundary, "Boundary conditions"); materials.clear(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModel::~SolidMechanicsModel() { AKANTU_DEBUG_IN(); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { delete *mat_it; } materials.clear(); + delete integrator; + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initVectors() { AKANTU_DEBUG_IN(); UInt nb_nodes = fem->getMesh().getNbNodes(); std::stringstream sstr_disp; sstr_disp << id << ":displacement"; std::stringstream sstr_mass; sstr_mass << id << ":mass"; std::stringstream sstr_velo; sstr_velo << id << ":velocity"; std::stringstream sstr_acce; sstr_acce << id << ":acceleration"; std::stringstream sstr_forc; sstr_forc << id << ":force"; std::stringstream sstr_resi; sstr_resi << id << ":residual"; std::stringstream sstr_boun; sstr_boun << id << ":boundary"; displacement = &(alloc(sstr_disp.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); mass = &(alloc(sstr_mass.str(), nb_nodes, 1)); // \todo see if it must not be spatial_dimension velocity = &(alloc(sstr_velo.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); acceleration = &(alloc(sstr_acce.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); force = &(alloc(sstr_forc.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); residual = &(alloc(sstr_resi.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); boundary = &(alloc(sstr_boun.str(), nb_nodes, spatial_dimension, false)); + std::stringstream sstr_curp; sstr_curp << id << ":current_position_tmp"; + current_position = &(alloc(sstr_curp.str(), 0, spatial_dimension, REAL_INIT_VALUE)); + const Mesh::ConnectivityTypeList & type_list = fem->getMesh().getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { UInt nb_element = fem->getMesh().getNbElement(*it); std::stringstream sstr_elma; sstr_elma << id << ":element_material:" << *it; element_material[*it] = &(alloc(sstr_elma.str(), nb_element, 1, 0)); } const Mesh::ConnectivityTypeList & ghost_type_list = fem->getMesh().getConnectivityTypeList(_ghost); for(it = ghost_type_list.begin(); it != ghost_type_list.end(); ++it) { UInt nb_element = fem->getMesh().getNbGhostElement(*it); std::stringstream sstr_elma; sstr_elma << id << ":ghost_element_material:" << *it; ghost_element_material[*it] = &(alloc(sstr_elma.str(), nb_element, 1, 0)); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::readMaterials(const std::string & filename) { std::ifstream infile; infile.open(filename.c_str()); std::string line; UInt current_line = 0, mat_count = materials.size(); if(!infile.good()) { AKANTU_DEBUG_ERROR("Canot open file " << filename); } while(infile.good()) { std::getline(infile, line); current_line++; /// remove comments size_t pos = line.find("#"); line = line.substr(0, pos); if(line.empty()) continue; std::stringstream sstr(line); std::string keyword; std::string value; sstr >> keyword; to_lower(keyword); if(keyword == "material") { std::string type; sstr >> type; to_lower(type); std::string obracket; sstr >> obracket; if(obracket != "[") AKANTU_DEBUG_ERROR("Malformed material file : missing [ at line " << current_line); std::stringstream sstr_mat; sstr_mat << id << ":" << mat_count++ << ":" << type; Material * material; MaterialID mat_id = sstr_mat.str(); if(type == "elastic") material = new MaterialElastic(*this, mat_id); else AKANTU_DEBUG_ERROR("Malformed material file : unknown material type " << type << " at line " << current_line); /// read the material properties std::getline(infile, line); size_t pos = line.find("#"); line = line.substr(0, pos); trim(line); current_line++; while(line[0] != ']') { pos = line.find("="); if(pos == std::string::npos) AKANTU_DEBUG_ERROR("Malformed material file : line must be \"key = value\" at line" << current_line); keyword = line.substr(0, pos); trim(keyword); value = line.substr(pos + 1); trim(value); try { material->setParam(keyword, value, mat_id); } catch (Exception ex) { AKANTU_DEBUG_ERROR("Malformed material file : error in setParam \"" << ex.info() << "\" at line " << current_line); } std::getline(infile, line); size_t pos = line.find("#"); line = line.substr(0, pos); trim(line); } materials.push_back(material); } } infile.close(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initMaterials() { Material ** mat_val = &(materials.at(0)); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { /// init internals properties (*mat_it)->initMaterial(); } /// 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); } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initModel() { fem->initShapeFunctions(_not_ghost); fem->initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMass() { AKANTU_DEBUG_IN(); UInt nb_nodes = fem->getMesh().getNbNodes(); memset(mass->values, 0, nb_nodes*sizeof(Real)); assembleMass(_not_ghost); assembleMass(_ghost); - /// @todo synchronize mass for the nodes of ghost elements - synchronize(_gst_smm_mass); - /// for not connected nodes put mass to one in order to avoid /// wrong range in paraview Real * mass_values = mass->values; for (UInt i = 0; i < nb_nodes; ++i) { if (!mass_values[i] || isnan(mass_values[i])) mass_values[i] = 1; } - + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMass(GhostType ghost_type) { AKANTU_DEBUG_IN(); Material ** mat_val = &(materials.at(0)); const Mesh::ConnectivityTypeList & type_list = fem->getMesh().getConnectivityTypeList(ghost_type); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); UInt nb_quadrature_points = FEM::getNbQuadraturePoints(*it); UInt shape_size = FEM::getShapeSize(*it); UInt nb_element; const Vector * shapes; UInt * elem_mat_val; if(ghost_type == _not_ghost) { nb_element = fem->getMesh().getNbElement(*it); shapes = &(fem->getShapes(*it)); elem_mat_val = element_material[*it]->values; } else { nb_element = fem->getMesh().getNbGhostElement(*it); shapes = &(fem->getGhostShapes(*it)); elem_mat_val = ghost_element_material[*it]->values; } if (nb_element == 0) continue; Vector * rho_phi_i = new Vector(nb_element, shape_size, "rho_x_shapes"); Real * rho_phi_i_val = rho_phi_i->values; Real * shapes_val = shapes->values; /// compute rho * \phi_i for each nodes of each element for (UInt el = 0; el < nb_element; ++el) { Real rho = mat_val[elem_mat_val[el]]->getRho(); for (UInt n = 0; n < shape_size; ++n) { *rho_phi_i_val++ = rho * *shapes_val++; } } Vector * int_rho_phi_i = new Vector(nb_element, shape_size / nb_quadrature_points, "inte_rho_x_shapes"); fem->integrate(*rho_phi_i, *int_rho_phi_i, nb_nodes_per_element, *it, ghost_type); delete rho_phi_i; fem->assembleVector(*int_rho_phi_i, *mass, 1, *it, ghost_type); delete int_rho_phi_i; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateResidual() { AKANTU_DEBUG_IN(); - /// @todo start synchronization - UInt nb_nodes = fem->getMesh().getNbNodes(); - Vector * current_position = new Vector(nb_nodes, spatial_dimension, NAN, "position"); + current_position->resize(nb_nodes); + //Vector * current_position = new Vector(nb_nodes, spatial_dimension, NAN, "position"); Real * current_position_val = current_position->values; Real * position_val = fem->getMesh().getNodes().values; Real * displacement_val = displacement->values; /// compute current_position = initial_position + displacement memcpy(current_position_val, position_val, nb_nodes*spatial_dimension*sizeof(Real)); for (UInt n = 0; n < nb_nodes*spatial_dimension; ++n) { *current_position_val++ += *displacement_val++; } + /// start synchronization + asynchronousSynchronize(_gst_smm_residual); + /// copy the forces in residual for boundary conditions memcpy(residual->values, force->values, nb_nodes*spatial_dimension*sizeof(Real)); /// call update residual on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->updateResidual(*current_position, _not_ghost); } - - /// @todo finalize communications + /// finalize communications + waitEndSynchronize(_gst_smm_residual); /// call update residual on each ghost elements for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->updateResidual(*current_position, _ghost); } - delete current_position; + // current_position; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateAcceleration() { AKANTU_DEBUG_IN(); UInt nb_nodes = acceleration->getSize(); + UInt nb_degre_of_freedom = acceleration->getNbComponent(); Real * mass_val = mass->values; Real * residual_val = residual->values; Real * accel_val = acceleration->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes; ++n) { for (UInt d = 0; d < nb_degre_of_freedom; d++) { if(!(*boundary_val)) { *accel_val = f_m2a * *residual_val / *mass_val; } residual_val++; accel_val++; boundary_val++; } mass_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitPred() { AKANTU_DEBUG_IN(); integrator->integrationSchemePred(time_step, *displacement, *velocity, *acceleration, *boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitCorr() { AKANTU_DEBUG_IN(); integrator->integrationSchemeCorr(time_step, *displacement, *velocity, *acceleration, *boundary); AKANTU_DEBUG_OUT(); } +/* -------------------------------------------------------------------------- */ +void SolidMechanicsModel::synchronizeBoundaries() { + AKANTU_DEBUG_IN(); + synchronize(_gst_smm_boundary); + AKANTU_DEBUG_OUT(); +}; + + /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep() { AKANTU_DEBUG_IN(); Material ** mat_val = &(materials.at(0)); Real min_dt = std::numeric_limits::max(); Real * coord = fem->getMesh().getNodes().values; Real * disp_val = displacement->values; const Mesh::ConnectivityTypeList & type_list = fem->getMesh().getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(fem->getMesh().getSpatialDimension(*it) != spatial_dimension) continue; UInt nb_nodes_per_element = fem->getMesh().getNbNodesPerElement(*it); UInt nb_element = fem->getMesh().getNbElement(*it); UInt * conn = fem->getMesh().getConnectivity(*it).values; UInt * elem_mat_val = element_material[*it]->values; Real * u = new Real[nb_nodes_per_element*spatial_dimension]; for (UInt el = 0; el < nb_element; ++el) { UInt el_offset = el * nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt offset_conn = conn[el_offset + n] * spatial_dimension; memcpy(u + n * spatial_dimension, coord + offset_conn, spatial_dimension * sizeof(Real)); for (UInt i = 0; i < spatial_dimension; ++i) { u[n * spatial_dimension + i] += disp_val[offset_conn + i]; } } Real el_size = fem->getElementInradius(u, *it); Real el_dt = mat_val[elem_mat_val[el]]->getStableTimeStep(el_size); min_dt = min_dt > el_dt ? el_dt : min_dt; } delete [] u; } - /// @todo reduction min over all processors + /// reduction min over all processors + allReduce(&min_dt, _so_min); + AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setPotentialEnergyFlagOn() { AKANTU_DEBUG_IN(); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->setPotentialEnergyFlagOn(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setPotentialEnergyFlagOff() { AKANTU_DEBUG_IN(); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->setPotentialEnergyFlagOff(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getPotentialEnergy() { AKANTU_DEBUG_IN(); Real epot = 0.; - /// fill the element filters of the materials using the element_material arrays - const Mesh::ConnectivityTypeList & type_list = fem->getMesh().getConnectivityTypeList(); - Mesh::ConnectivityTypeList::const_iterator it; - for(it = type_list.begin(); it != type_list.end(); ++it) { - if(fem->getMesh().getSpatialDimension(*it) != spatial_dimension) continue; - - std::vector::iterator mat_it; - for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { - epot += fem->integrate((*mat_it)->getPotentialEnergy(*it), - *it); - } + /// call update residual on each local elements + std::vector::iterator mat_it; + for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { + epot += (*mat_it)->getPotentialEnergy(); } - /// @todo reduction sum over all processors + /// reduction sum over all processors + allReduce(&epot, _so_sum); AKANTU_DEBUG_OUT(); return epot; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy() { AKANTU_DEBUG_IN(); - UInt nb_nodes = velocity->getSize(); - UInt nb_degre_of_freedom = velocity->getNbComponent(); - - Real * mass_val = mass->values; - Real * vel_val = velocity->values; Real ekin = 0.; + UInt nb_nodes = fem->getMesh().getNbNodes(); + Vector * v_square = new Vector(nb_nodes, 1, "v_square"); + + Real * vel_val = velocity->values; + Real * v_s_val = v_square->values; + for (UInt n = 0; n < nb_nodes; ++n) { - Real norm_vel = 0.; - for (UInt d = 0; d < nb_degre_of_freedom; d++) { - norm_vel += *vel_val * *vel_val; + *v_s_val = 0; + for (UInt s = 0; s < spatial_dimension; ++s) { + *v_s_val += *vel_val * *vel_val; vel_val++; } - ekin += *mass_val * norm_vel; + v_s_val++; + } - mass_val++; + Material ** mat_val = &(materials.at(0)); + + const Mesh:: ConnectivityTypeList & type_list = fem->getMesh().getConnectivityTypeList(); + Mesh::ConnectivityTypeList::const_iterator it; + for(it = type_list.begin(); it != type_list.end(); ++it) { + if(fem->getMesh().getSpatialDimension(*it) != spatial_dimension) continue; + + UInt nb_quadrature_points = FEM::getNbQuadraturePoints(*it); + UInt nb_element = fem->getMesh().getNbElement(*it); + + Vector * v_square_el = new Vector(nb_element, nb_quadrature_points, "v_square per element"); + + fem->interpolateOnQuadraturePoints(*v_square, *v_square_el, 1, *it); + + Real * v_square_el_val = v_square_el->values; + UInt * elem_mat_val = element_material[*it]->values; + + for (UInt el = 0; el < nb_element; ++el) { + Real rho = mat_val[elem_mat_val[el]]->getRho(); + for (UInt q = 0; q < nb_quadrature_points; ++q) { + *v_square_el_val *= rho; + v_square_el_val++; + } + } + + ekin += fem->integrate(*v_square_el, *it); + + delete v_square_el; } + delete v_square; - /// @todo reduction sum over all processors + /// reduction sum over all processors + allReduce(&ekin, _so_sum); AKANTU_DEBUG_OUT(); return ekin * .5; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "Solid Mechanics Model [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + spatial dimension : " << spatial_dimension << std::endl; stream << space << " + fem [" << std::endl; fem->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + nodals information [" << std::endl; displacement->printself(stream, indent + 2); mass ->printself(stream, indent + 2); velocity ->printself(stream, indent + 2); acceleration->printself(stream, indent + 2); force ->printself(stream, indent + 2); residual ->printself(stream, indent + 2); boundary ->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + connectivity type information [" << std::endl; const Mesh::ConnectivityTypeList & type_list = fem->getMesh().getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { stream << space << AKANTU_INDENT << AKANTU_INDENT << " + " << *it <<" [" << std::endl; element_material[*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/model/solid_mechanics/solid_mechanics_model.hh b/model/solid_mechanics/solid_mechanics_model.hh index ab5cbd5e5..08b70dce0 100644 --- a/model/solid_mechanics/solid_mechanics_model.hh +++ b/model/solid_mechanics/solid_mechanics_model.hh @@ -1,206 +1,212 @@ /** * @file solid_mechanics_model.hh * @author Nicolas Richart * @date Thu Jul 22 11:51:06 2010 * * @brief Model of Solid Mechanics * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_MECHANICS_MODEL_HH__ #define __AKANTU_SOLID_MECHANICS_MODEL_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "model.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ namespace akantu { // class Material; class IntegrationScheme2ndOrder; } __BEGIN_AKANTU__ class SolidMechanicsModel : public Model { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: SolidMechanicsModel(UInt spatial_dimension, const ModelID & id = "solid_mechanics_model", const MemoryID & memory_id = 0); SolidMechanicsModel(Mesh & mesh, UInt spatial_dimension = 0, const ModelID & id = "solid_mechanics_model", const MemoryID & memory_id = 0); virtual ~SolidMechanicsModel(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// allocate all vectors void initVectors(); /// read the material files to instantiate all the materials void readMaterials(const std::string & filename); /// initialize all internal arrays for materials void initMaterials(); /// initialize the model void initModel(); /// assemble the lumped mass matrix void assembleMass(); /// assemble the residual for the explicit scheme void updateResidual(); /// compute the acceleration from the residual void updateAcceleration(); /// explicit integration predictor void explicitPred(); /// explicit integration corrector void explicitCorr(); /// compute boundary forces from quadrature point force values void computeForcesFromQuadraturePointForceValues(); + /// synchronize the ghost element boundaries values + void synchronizeBoundaries(); + /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; private: /// assemble the lumped mass matrix for local and ghost elements void assembleMass(GhostType ghost_type); /* ------------------------------------------------------------------------ */ /* Ghost Synchronizer inherited members */ /* ------------------------------------------------------------------------ */ public: inline virtual UInt getNbDataToPack(const Element & element, GhostSynchronizationTag tag) const; inline virtual UInt getNbDataToUnpack(const Element & element, GhostSynchronizationTag tag) const; inline virtual void packData(Real ** buffer, const Element & element, GhostSynchronizationTag tag) const; inline virtual void unpackData(Real ** buffer, const Element & element, GhostSynchronizationTag tag) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(TimeStep, time_step, Real); AKANTU_SET_MACRO(TimeStep, time_step, Real); AKANTU_GET_MACRO(F_M2A, f_m2a, Real); AKANTU_SET_MACRO(F_M2A, f_m2a, Real); AKANTU_GET_MACRO(Displacement, *displacement, Vector &); AKANTU_GET_MACRO(Mass, *mass, Vector &); AKANTU_GET_MACRO(Velocity, *velocity, Vector &); AKANTU_GET_MACRO(Acceleration, *acceleration, Vector &); AKANTU_GET_MACRO(Force, *force, Vector &); AKANTU_GET_MACRO(Residual, *residual, Vector &); AKANTU_GET_MACRO(Boundary, *boundary, Vector &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementMaterial, element_material, Vector &); inline Material & getMaterial(UInt mat_index); /// compute the stable time step Real getStableTimeStep(); void setPotentialEnergyFlagOn(); void setPotentialEnergyFlagOff(); Real getPotentialEnergy(); Real getKineticEnergy(); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// time step Real time_step; /// conversion coefficient form force/mass to acceleration Real f_m2a; /// displacements array Vector * displacement; /// lumped mass array Vector * mass; /// velocities array Vector * velocity; /// accelerations array Vector * acceleration; /// forces array Vector * force; /// residuals array Vector * residual; /// boundaries array Vector * boundary; + /// array of current position used during update residual + Vector * current_position; + /// materials of all elements ByElementTypeUInt element_material; /// materials of all ghost elements ByElementTypeUInt ghost_element_material; /// list of used materials std::vector materials; /// integration scheme of second order used IntegrationScheme2ndOrder * integrator; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model_inline_impl.cc" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const SolidMechanicsModel & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_SOLID_MECHANICS_MODEL_HH__ */ diff --git a/model/solid_mechanics/solid_mechanics_model_inline_impl.cc b/model/solid_mechanics/solid_mechanics_model_inline_impl.cc index 04ec2f646..d8889d23a 100644 --- a/model/solid_mechanics/solid_mechanics_model_inline_impl.cc +++ b/model/solid_mechanics/solid_mechanics_model_inline_impl.cc @@ -1,185 +1,225 @@ /** * @file solid_mechanics_model_inline_impl.cc * @author Nicolas Richart * @date Thu Jul 29 12:07:04 2010 * * @brief Implementation of the inline functions of the SolidMechanicsModel class * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ inline Material & SolidMechanicsModel::getMaterial(UInt mat_index) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(mat_index < materials.size(), "The model " << id << " has no material no "<< mat_index); AKANTU_DEBUG_OUT(); return *materials[mat_index]; } /* -------------------------------------------------------------------------- */ inline UInt SolidMechanicsModel::getNbDataToPack(const Element & element, GhostSynchronizationTag tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); #ifdef AKANTU_DEBUG - size += spatial_dimension; /// position of the barycenter of the P1 element (only for check) + size += spatial_dimension; /// position of the barycenter of the element (only for check) #endif switch(tag) { case _gst_smm_mass: { size += nb_nodes_per_element; // mass vector break; } case _gst_smm_residual: { + size += nb_nodes_per_element * spatial_dimension; // displacement + UInt mat = element_material[element.type]->values[element.element]; size += materials[mat]->getNbDataToPack(element, tag); break; } + case _gst_smm_boundary: { + size += nb_nodes_per_element * spatial_dimension * 3; // force, displacement, boundary + break; + } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; }; /* -------------------------------------------------------------------------- */ inline UInt SolidMechanicsModel::getNbDataToUnpack(const Element & element, GhostSynchronizationTag tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); #ifdef AKANTU_DEBUG - size += spatial_dimension; /// position of the barycenter of the P1 element (only for check) + size += spatial_dimension; /// position of the barycenter of the element (only for check) #endif switch(tag) { case _gst_smm_mass: { size += nb_nodes_per_element; // mass vector break; } case _gst_smm_residual: { + size += nb_nodes_per_element * spatial_dimension; // displacement + UInt mat = ghost_element_material[element.type]->values[element.element]; size += materials[mat]->getNbDataToPack(element, tag); break; } + case _gst_smm_boundary: { + size += nb_nodes_per_element * spatial_dimension * 3; // force, displacement, boundary + break; + } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); - return 0; + return size; }; /* -------------------------------------------------------------------------- */ inline void SolidMechanicsModel::packData(Real ** buffer, const Element & element, GhostSynchronizationTag tag) const { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); UInt el_offset = element.element * nb_nodes_per_element; - -#ifdef AKANTU_DEBUG - Real * coord = fem->getMesh().getNodes().values; UInt * conn = fem->getMesh().getConnectivity(element.type).values; - UInt spatial_dimension = Mesh::getSpatialDimension(element.type); - memset(*buffer, 0, spatial_dimension * sizeof(Real)); - for (UInt n = 0; n < nb_nodes_per_element; ++n) { - UInt offset_conn = conn[el_offset + n] * spatial_dimension; - for (UInt i = 0; i < spatial_dimension; ++i) { - *buffer[i] += coord[offset_conn + i] / (Real) nb_nodes_per_element; - } - } - *buffer += spatial_dimension; +#ifdef AKANTU_DEBUG + fem->getMesh().getBarycenter(element.element, element.type, *buffer); + (*buffer) += spatial_dimension; #endif - switch(tag) { case _gst_smm_mass: { for (UInt n = 0; n < nb_nodes_per_element; ++n) { - UInt offset_conn = conn[el_offset + n] * spatial_dimension; - **buffer = mass->values[offset_conn]; - *buffer++; + UInt offset_conn = conn[el_offset + n]; + (*buffer)[n] = mass->values[offset_conn]; } + *buffer += nb_nodes_per_element; break; } case _gst_smm_residual: { - UInt mat = ghost_element_material[element.type]->values[element.element]; + for (UInt n = 0; n < nb_nodes_per_element; ++n) { + UInt offset_conn = conn[el_offset + n] * spatial_dimension; + memcpy(*buffer, current_position->values + offset_conn, spatial_dimension * sizeof(Real)); + *buffer += spatial_dimension; + } + + UInt mat = element_material[element.type]->values[element.element]; materials[mat]->packData(buffer, element, tag); break; } + case _gst_smm_boundary: { + for (UInt n = 0; n < nb_nodes_per_element; ++n) { + UInt offset_conn = conn[el_offset + n] * spatial_dimension; + + memcpy(*buffer, force->values + offset_conn, spatial_dimension * sizeof(Real)); + *buffer += spatial_dimension; + + memcpy(*buffer, velocity->values + offset_conn, spatial_dimension * sizeof(Real)); + *buffer += spatial_dimension; + + for (UInt i = 0; i < spatial_dimension; ++i) { + (*buffer)[i] = boundary->values[offset_conn + i] ? 1.0 : -1.0; + } + *buffer += spatial_dimension; + } + break; + } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); }; /* -------------------------------------------------------------------------- */ inline void SolidMechanicsModel::unpackData(Real ** buffer, const Element & element, GhostSynchronizationTag tag) const { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); UInt el_offset = element.element * nb_nodes_per_element; + UInt * conn = fem->getMesh().getGhostConnectivity(element.type).values; #ifdef AKANTU_DEBUG - Real * coord = fem->getMesh().getNodes().values; - UInt * conn = fem->getMesh().getConnectivity(element.type).values; - UInt spatial_dimension = Mesh::getSpatialDimension(element.type); - Real barycenter[spatial_dimension]; - memset(barycenter, 0, spatial_dimension * sizeof(Real)); - for (UInt n = 0; n < nb_nodes_per_element; ++n) { - UInt offset_conn = conn[el_offset + n] * spatial_dimension; - for (UInt i = 0; i < spatial_dimension; ++i) { - barycenter[i] += coord[offset_conn + i] / (Real) nb_nodes_per_element; - } - } + fem->getMesh().getBarycenter(element.element, element.type, barycenter, _ghost); Real tolerance = 1e-15; - for (UInt i = 0; i < spatial_dimension; ++i) { - if(abs(barycenter[i] - *buffer[i]) <= tolerance) + if(!(fabs(barycenter[i] - (*buffer)[i]) <= tolerance)) AKANTU_DEBUG_ERROR("Unpacking an unknown value for the element : " - << element); + << element + << "(barycenter[" << i << "] = " << barycenter[i] + << " and (*buffer)[" << i << "] = " << (*buffer)[i] << ")"); } *buffer += spatial_dimension; #endif switch(tag) { case _gst_smm_mass: { for (UInt n = 0; n < nb_nodes_per_element; ++n) { - UInt offset_conn = conn[el_offset + n] * spatial_dimension; - mass->values[offset_conn] = **buffer; - *buffer++; + UInt offset_conn = conn[el_offset + n]; + mass->values[offset_conn] = (*buffer)[n]; } + *buffer += nb_nodes_per_element; break; } case _gst_smm_residual: { + for (UInt n = 0; n < nb_nodes_per_element; ++n) { + UInt offset_conn = conn[el_offset + n] * spatial_dimension; + memcpy(current_position->values + offset_conn, *buffer, spatial_dimension * sizeof(Real)); + *buffer += spatial_dimension; + } + UInt mat = ghost_element_material[element.type]->values[element.element]; materials[mat]->unpackData(buffer, element, tag); break; } + case _gst_smm_boundary: { + for (UInt n = 0; n < nb_nodes_per_element; ++n) { + UInt offset_conn = conn[el_offset + n] * spatial_dimension; + + memcpy(force->values + offset_conn, *buffer, spatial_dimension * sizeof(Real)); + *buffer += spatial_dimension; + + memcpy(velocity->values + offset_conn, *buffer, spatial_dimension * sizeof(Real)); + *buffer += spatial_dimension; + + for (UInt i = 0; i < spatial_dimension; ++i) { + boundary->values[offset_conn + i] = ((*buffer)[i] > 0) ? true : false; + } + *buffer += spatial_dimension; + } + break; + } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); }; diff --git a/synchronizer/communicator.cc b/synchronizer/communicator.cc index 6a34c6eca..9a9ca8968 100644 --- a/synchronizer/communicator.cc +++ b/synchronizer/communicator.cc @@ -1,579 +1,580 @@ /** * @file communicator.cc * @author Nicolas Richart * @date Thu Aug 26 16:36:21 2010 * * @brief implementation of a communicator using a static_communicator for real * send/receive * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "ghost_synchronizer.hh" #include "communicator.hh" #include "static_communicator.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ Communicator::Communicator(SynchronizerID id, MemoryID memory_id) : Synchronizer(id, memory_id), static_communicator(StaticCommunicator::getStaticCommunicator()) { AKANTU_DEBUG_IN(); for (UInt t = _not_defined; t < _max_element_type; ++t) { element_to_send_offset [t] = NULL; element_to_send [t] = NULL; element_to_receive_offset[t] = NULL; element_to_receive [t] = NULL; } - UInt nb_proc = static_communicator->getNbProc(); + nb_proc = static_communicator->getNbProc(); + rank = static_communicator->whoAmI(); + send_buffer = new Vector[nb_proc]; recv_buffer = new Vector[nb_proc]; send_element = new std::vector[nb_proc]; recv_element = new std::vector[nb_proc]; - size_to_send.clear(); size_to_receive.clear(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Communicator::~Communicator() { AKANTU_DEBUG_IN(); - UInt psize = static_communicator->getNbProc(); - for (UInt p = 0; p < psize; ++p) { + for (UInt p = 0; p < nb_proc; ++p) { send_buffer[p].resize(0); recv_buffer[p].resize(0); send_element[p].clear(); recv_element[p].clear(); } delete [] send_buffer; delete [] recv_buffer; delete [] send_element; delete [] recv_element; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Communicator * Communicator::createCommunicatorDistributeMesh(Mesh & mesh, const MeshPartition * partition, UInt root, SynchronizerID id, MemoryID memory_id) { AKANTU_DEBUG_IN(); StaticCommunicator * comm = StaticCommunicator::getStaticCommunicator(); - UInt nb_proc = comm->getNbProc(); UInt my_rank = comm->whoAmI(); UInt * local_connectivity; UInt * local_partitions; Vector old_nodes; Vector * nodes = mesh.getNodesPointer(); UInt spatial_dimension = nodes->getNbComponent(); Communicator * communicator = new Communicator(); if(nb_proc == 1) return communicator; /* ------------------------------------------------------------------------ */ /* Local (rank == root) */ /* ------------------------------------------------------------------------ */ if(my_rank == root) { AKANTU_DEBUG_ASSERT(partition->getNbPartition() == nb_proc, "The number of partition does not match the number of processors"); /** * connectivity and communications scheme construction */ 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; if(Mesh::getSpatialDimension(type) != mesh.getSpatialDimension()) continue; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_element = mesh.getNbElement(*it); UInt nb_local_element[nb_proc]; UInt nb_ghost_element[nb_proc]; UInt nb_element_to_send[nb_proc]; memset(nb_local_element, 0, nb_proc*sizeof(UInt)); memset(nb_ghost_element, 0, nb_proc*sizeof(UInt)); memset(nb_element_to_send, 0, nb_proc*sizeof(UInt)); UInt * partition_num = partition->getPartition(type).values; UInt * ghost_partition = partition->getGhostPartition(type).values; UInt * ghost_partition_offset = partition->getGhostPartitionOffset(type).values; /// constricting the reordering structures for (UInt el = 0; el < nb_element; ++el) { nb_local_element[partition_num[el]]++; for (UInt part = ghost_partition_offset[el]; part < ghost_partition_offset[el + 1]; ++part) { nb_ghost_element[ghost_partition[part]]++; } nb_element_to_send[partition_num[el]] += ghost_partition_offset[el + 1] - ghost_partition_offset[el] + 1; } /// allocating buffers UInt * buffers[nb_proc]; UInt * buffers_tmp[nb_proc]; for (UInt p = 0; p < nb_proc; ++p) { UInt size = nb_nodes_per_element * (nb_local_element[p] + nb_ghost_element[p]); buffers[p] = new UInt[size]; buffers_tmp[p] = buffers[p]; } /// copying the local connectivity UInt * conn_val = mesh.getConnectivity(type).values; for (UInt el = 0; el < nb_element; ++el) { memcpy(buffers_tmp[partition_num[el]], conn_val + el * nb_nodes_per_element, nb_nodes_per_element * sizeof(UInt)); buffers_tmp[partition_num[el]] += nb_nodes_per_element; } /// copying the connectivity of ghost element for (UInt el = 0; el < nb_element; ++el) { for (UInt part = ghost_partition_offset[el]; part < ghost_partition_offset[el + 1]; ++part) { UInt proc = ghost_partition[part]; memcpy(buffers_tmp[proc], conn_val + el * nb_nodes_per_element, nb_nodes_per_element * sizeof(UInt)); buffers_tmp[proc] += nb_nodes_per_element; } } /// send all connectivity and ghost information to all processors std::vector requests; for (UInt p = 0; p < nb_proc; ++p) { if(p != root) { UInt size[4]; size[0] = (UInt) type; size[1] = nb_local_element[p]; size[2] = nb_ghost_element[p]; size[3] = nb_element_to_send[p]; comm->send(size, 4, p, 0); AKANTU_DEBUG_INFO("Sending connectivities to proc " << p); requests.push_back(comm->asyncSend(buffers[p], nb_nodes_per_element * (nb_local_element[p] + nb_ghost_element[p]), p, 1)); } else { local_connectivity = buffers[p]; } } /// create the renumbered connectivity AKANTU_DEBUG_INFO("Renumbering local connectivities"); MeshUtils::renumberMeshNodes(mesh, local_connectivity, nb_local_element[root], nb_ghost_element[root], type, old_nodes); comm->waitAll(requests); comm->freeCommunicationRequest(requests); requests.clear(); for (UInt p = 0; p < nb_proc; ++p) { delete [] buffers[p]; } for (UInt p = 0; p < nb_proc; ++p) { buffers[p] = new UInt[nb_ghost_element[p] + nb_element_to_send[p]]; buffers_tmp[p] = buffers[p]; } /// splitting the partition information to send them to processors UInt count_by_proc[nb_proc]; memset(count_by_proc, 0, nb_proc*sizeof(UInt)); for (UInt el = 0; el < nb_element; ++el) { *(buffers_tmp[partition_num[el]]++) = ghost_partition_offset[el + 1] - ghost_partition_offset[el]; for (UInt part = ghost_partition_offset[el], i = 0; part < ghost_partition_offset[el + 1]; ++part, ++i) { *(buffers_tmp[partition_num[el]]++) = ghost_partition[part]; } } for (UInt el = 0; el < nb_element; ++el) { for (UInt part = ghost_partition_offset[el], i = 0; part < ghost_partition_offset[el + 1]; ++part, ++i) { *(buffers_tmp[ghost_partition[part]]++) = partition_num[el]; } } /// last data to compute the communication scheme for (UInt p = 0; p < nb_proc; ++p) { if(p != root) { AKANTU_DEBUG_INFO("Sending partition informations to proc " << p); requests.push_back(comm->asyncSend(buffers[p], nb_element_to_send[p] + nb_ghost_element[p], p, 2)); } else { local_partitions = buffers[p]; } } AKANTU_DEBUG_INFO("Creating communications scheme"); communicator->fillCommunicationScheme(local_partitions, nb_local_element[root], nb_ghost_element[root], nb_element_to_send[root], type); comm->waitAll(requests); comm->freeCommunicationRequest(requests); requests.clear(); for (UInt p = 0; p < nb_proc; ++p) { delete [] buffers[p]; } } for (UInt p = 0; p < nb_proc; ++p) { if(p != root) { UInt size[4]; size[0] = (UInt) _not_defined; size[1] = 0; size[2] = 0; size[3] = 0; comm->send(size, 4, p, 0); } } /** * Nodes coordinate construction and synchronization */ /// get the list of nodes to send and send them Real * local_nodes; for (UInt p = 0; p < nb_proc; ++p) { UInt nb_nodes; UInt * buffer; if(p != root) { AKANTU_DEBUG_INFO("Receiving list of nodes from proc " << p); comm->receive(&nb_nodes, 1, p, 0); buffer = new UInt[nb_nodes]; comm->receive(buffer, nb_nodes, p, 3); } else { nb_nodes = old_nodes.getSize(); buffer = old_nodes.values; } /// get the coordinates for the selected nodes Real * nodes_to_send = new Real[nb_nodes * spatial_dimension]; Real * nodes_to_send_tmp = nodes_to_send; for (UInt n = 0; n < nb_nodes; ++n) { memcpy(nodes_to_send_tmp, nodes->values + spatial_dimension * buffer[n], spatial_dimension * sizeof(Real)); nodes_to_send_tmp += spatial_dimension; } if(p != root) { /// send them for distant processors delete [] buffer; AKANTU_DEBUG_INFO("Sending coordinates to proc " << p); comm->send(nodes_to_send, nb_nodes * spatial_dimension, p, 4); delete [] nodes_to_send; } else { /// save them for local processor local_nodes = nodes_to_send; } } /// construct the local nodes coordinates UInt nb_nodes = old_nodes.getSize(); nodes->resize(nb_nodes); memcpy(nodes->values, local_nodes, nb_nodes * spatial_dimension * sizeof(Real)); delete [] local_nodes; /* ---------------------------------------------------------------------- */ /* Distant (rank != root) */ /* ---------------------------------------------------------------------- */ } else { /** * connectivity and communications scheme construction on distant processors */ ElementType type = _not_defined; do { UInt size[4]; comm->receive(size, 4, root, 0); type = (ElementType) size[0]; UInt nb_local_element = size[1]; UInt nb_ghost_element = size[2]; UInt nb_element_to_send = size[3]; if(type != _not_defined) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); local_connectivity = new UInt[(nb_local_element + nb_ghost_element) * nb_nodes_per_element]; AKANTU_DEBUG_INFO("Receiving connectivities from proc " << root); comm->receive(local_connectivity, nb_nodes_per_element * (nb_local_element + nb_ghost_element), root, 1); AKANTU_DEBUG_INFO("Renumbering local connectivities"); MeshUtils::renumberMeshNodes(mesh, local_connectivity, nb_local_element, nb_ghost_element, type, old_nodes); delete [] local_connectivity; local_partitions = new UInt[nb_element_to_send + nb_ghost_element * 2]; AKANTU_DEBUG_INFO("Receiving partition informations from proc " << root); comm->receive(local_partitions, nb_element_to_send + nb_ghost_element * 2, root, 2); AKANTU_DEBUG_INFO("Creating communications scheme"); communicator->fillCommunicationScheme(local_partitions, nb_local_element, nb_ghost_element, nb_element_to_send, type); delete [] local_partitions; } } while(type != _not_defined); /** * Nodes coordinate construction and synchronization on distant processors */ AKANTU_DEBUG_INFO("Sending list of nodes to proc " << root); UInt nb_nodes = old_nodes.getSize(); comm->send(&nb_nodes, 1, root, 0); comm->send(old_nodes.values, nb_nodes, root, 3); nodes->resize(nb_nodes); AKANTU_DEBUG_INFO("Receiving coordinates from proc " << root); comm->receive(nodes->values, nb_nodes * spatial_dimension, root, 4); } AKANTU_DEBUG_OUT(); return communicator; } /* -------------------------------------------------------------------------- */ void Communicator::fillCommunicationScheme(UInt * partition, UInt nb_local_element, UInt nb_ghost_element, UInt nb_element_to_send, ElementType type) { AKANTU_DEBUG_IN(); Element element; element.type = type; UInt * part = partition; part = partition; for (UInt lel = 0; lel < nb_local_element; ++lel) { UInt nb_send = *part; part++; element.element = lel; for (UInt p = 0; p < nb_send; ++p) { UInt proc = *part; part++; AKANTU_DEBUG(dblAccessory, "Must send : " << element << " to proc " << proc); (send_element[proc]).push_back(element); } } for (UInt gel = 0; gel < nb_ghost_element; ++gel) { UInt proc = *part; part++; element.element = gel; AKANTU_DEBUG(dblAccessory, "Must recv : " << element << " from proc " << proc); recv_element[proc].push_back(element); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Communicator::synchronize(GhostSynchronizationTag tag) { AKANTU_DEBUG_IN(); asynchronousSynchronize(tag); waitEndSynchronize(tag); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Communicator::asynchronousSynchronize(GhostSynchronizationTag tag) { AKANTU_DEBUG_IN(); - UInt psize = static_communicator->getNbProc(); - UInt prank = static_communicator->whoAmI(); - AKANTU_DEBUG_ASSERT(send_requests.size() == 0, "There must be some pending sending communications"); - for (UInt p = 0; p < psize; ++p) { + for (UInt p = 0; p < nb_proc; ++p) { UInt ssize = (size_to_send[tag]).values[p]; - if(p == prank || ssize == 0) continue; + if(p == rank || ssize == 0) continue; - AKANTU_DEBUG_INFO("Packing data for proc" << p); send_buffer[p].resize(ssize); Real * buffer = send_buffer[p].values; - Element * elements = &(send_element[p].at(0)); UInt nb_elements = send_element[p].size(); + AKANTU_DEBUG_INFO("Packing data for proc " << p + << " (" << ssize << "/" << nb_elements + <<" data to send/elements)"); for (UInt el = 0; el < nb_elements; ++el) { ghost_synchronizer->packData(&buffer, *elements, tag); elements++; } - + std::cerr << std::dec; AKANTU_DEBUG_INFO("Posting send to proc " << p); send_requests.push_back(static_communicator->asyncSend(send_buffer[p].values, ssize, p, (Int) tag)); } AKANTU_DEBUG_ASSERT(recv_requests.size() == 0, "There must be some pending receive communications"); - for (UInt p = 0; p < psize; ++p) { - UInt rsize = size_to_receive[tag].values[p]; - if(p == prank || rsize == 0) continue; + for (UInt p = 0; p < nb_proc; ++p) { + UInt rsize = (size_to_receive[tag]).values[p]; + if(p == rank || rsize == 0) continue; recv_buffer[p].resize(rsize); - AKANTU_DEBUG_INFO("Posting receive from proc " << p); + AKANTU_DEBUG_INFO("Posting receive from proc " << p << " (" << rsize << " data to receive)"); recv_requests.push_back(static_communicator->asyncReceive(recv_buffer[p].values, rsize, p, (Int) tag)); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Communicator::waitEndSynchronize(GhostSynchronizationTag tag) { AKANTU_DEBUG_IN(); std::vector req_not_finished; std::vector * req_not_finished_tmp = &req_not_finished; std::vector * recv_requests_tmp = &recv_requests; while(!recv_requests_tmp->empty()) { for (std::vector::iterator req_it = recv_requests_tmp->begin(); req_it != recv_requests_tmp->end() ; ++req_it) { CommunicationRequest * req = *req_it; if(static_communicator->testRequest(req)) { UInt proc = req->getSource(); AKANTU_DEBUG_INFO("Unpacking data coming from proc " << proc); Real * buffer = recv_buffer[proc].values; Element * elements = &recv_element[proc].at(0); UInt nb_elements = recv_element[proc].size(); for (UInt el = 0; el < nb_elements; ++el) { ghost_synchronizer->unpackData(&buffer, *elements, tag); elements++; } static_communicator->freeCommunicationRequest(req); } else { req_not_finished_tmp->push_back(req); } } std::vector * swap = req_not_finished_tmp; req_not_finished_tmp = recv_requests_tmp; recv_requests_tmp = swap; req_not_finished_tmp->clear(); } static_communicator->waitAll(send_requests); static_communicator->freeCommunicationRequest(send_requests); send_requests.clear(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -void Communicator::registerTag(GhostSynchronizationTag tag) { +void Communicator::allReduce(Real * values, UInt nb_values, const SynchronizerOperation & op) { AKANTU_DEBUG_IN(); + static_communicator->allReduce(values, nb_values, op); + AKANTU_DEBUG_OUT(); +} - UInt psize = static_communicator->getNbProc(); - UInt prank = static_communicator->whoAmI(); +/* -------------------------------------------------------------------------- */ +void Communicator::registerTag(GhostSynchronizationTag tag) { + AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(size_to_send.find(tag) == size_to_send.end(), "The GhostSynchronizationTag " << tag << "is already registered in " << id); - size_to_send [tag].resize(psize); - size_to_receive[tag].resize(psize); + size_to_send [tag].resize(nb_proc); + size_to_receive[tag].resize(nb_proc); - for (UInt p = 0; p < psize; ++p) { + for (UInt p = 0; p < nb_proc; ++p) { UInt ssend = 0; UInt sreceive = 0; - if(p != prank) { + if(p != rank) { for (std::vector::const_iterator sit = send_element[p].begin(); sit != send_element[p].end(); ++sit) { ssend += ghost_synchronizer->getNbDataToPack(*sit, tag); } for (std::vector::const_iterator rit = recv_element[p].begin(); rit != recv_element[p].end(); ++rit) { sreceive += ghost_synchronizer->getNbDataToUnpack(*rit, tag); } AKANTU_DEBUG_INFO("I have " << ssend << " data to send to " << p - << " and " << sreceive << " data to receive"); + << " and " << sreceive << " data to receive for tag " << tag); } size_to_send [tag].values[p] = ssend; size_to_receive[tag].values[p] = sreceive; } AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/synchronizer/communicator.hh b/synchronizer/communicator.hh index 2b10d8855..a490c92f0 100644 --- a/synchronizer/communicator.hh +++ b/synchronizer/communicator.hh @@ -1,121 +1,128 @@ /** * @file communicator.hh * @author Nicolas Richart * @date Thu Aug 19 15:28:35 2010 * * @brief wrapper to the static communicator * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COMMUNICATOR_HH__ #define __AKANTU_COMMUNICATOR_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_vector.hh" #include "static_communicator.hh" #include "synchronizer.hh" #include "mesh.hh" #include "mesh_partition.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ class Communicator : public Synchronizer { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Communicator(SynchronizerID id = "communicator", MemoryID memory_id = 0); virtual ~Communicator(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// get a mesh and a partition and create the local mesh and the associated communicator static Communicator * createCommunicatorDistributeMesh(Mesh & mesh, const MeshPartition * partition, UInt root = 0, SynchronizerID id = "communicator", MemoryID memory_id = 0); /* ------------------------------------------------------------------------ */ /* Inherited from Synchronizer */ /* ------------------------------------------------------------------------ */ /// synchronize ghosts void synchronize(GhostSynchronizationTag tag); /// asynchronous synchronization of ghosts void asynchronousSynchronize(GhostSynchronizationTag tag); /// wait end of asynchronous synchronization of ghosts void waitEndSynchronize(GhostSynchronizationTag tag); + /// do a all reduce operation + void allReduce(Real * values, UInt nb_values, const SynchronizerOperation & op); + + /* ------------------------------------------------------------------------ */ /// register a new communication void registerTag(GhostSynchronizationTag tag); protected: /// fill the communications array of a communicator based on a partition array void fillCommunicationScheme(UInt * partition, UInt nb_local_element, UInt nb_ghost_element, UInt nb_element_to_send, ElementType type); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// the static memory instance StaticCommunicator * static_communicator; /// size of data to send to each processor by communication tag std::map< GhostSynchronizationTag, Vector > size_to_send; /// size of data to receive form each processor by communication tag std::map< GhostSynchronizationTag, Vector > size_to_receive; Vector * send_buffer; Vector * recv_buffer; /// send requests std::vector send_requests; /// receive requests std::vector recv_requests; /// list of real element to send ordered by type then by receiver processors ByElementTypeUInt element_to_send_offset; ByElementTypeUInt element_to_send; std::vector * send_element; std::vector * recv_element; /// list of ghost element to receive ordered by type then by sender processors ByElementTypeUInt element_to_receive_offset; ByElementTypeUInt element_to_receive; + + UInt nb_proc; + UInt rank; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "communicator_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_COMMUNICATOR_HH__ */ diff --git a/synchronizer/ghost_synchronizer.cc b/synchronizer/ghost_synchronizer.cc index f5e4d154d..7113fd8b5 100644 --- a/synchronizer/ghost_synchronizer.cc +++ b/synchronizer/ghost_synchronizer.cc @@ -1,135 +1,150 @@ /** * @file ghost_synchronizer.cc * @author Nicolas Richart * @date Tue Aug 24 18:59:17 2010 * * @brief implementation of the ghost synchronizer * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #include "ghost_synchronizer.hh" #include "synchronizer.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ GhostSynchronizer::GhostSynchronizer() : nb_ghost_synchronization_tags(0) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ GhostSynchronizer::~GhostSynchronizer() { AKANTU_DEBUG_IN(); for (std::list::iterator it = synchronizers.begin(); it != synchronizers.end(); ++it) { delete (*it); } synchronizers.clear(); registered_synchronization.clear(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void GhostSynchronizer::synchronize(GhostSynchronizationTag tag) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(registered_synchronization.find(tag) != registered_synchronization.end(), "Tag " << tag << " is not registered."); AKANTU_DEBUG_INFO("Synchronizing the tag : " << registered_synchronization.find(tag)->second << " (" << tag <<")"); for (std::list::iterator it = synchronizers.begin(); it != synchronizers.end(); ++it) { (*it)->synchronize(tag); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void GhostSynchronizer::asynchronousSynchronize(GhostSynchronizationTag tag) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(registered_synchronization.find(tag) != registered_synchronization.end(), "Tag " << tag << " is not registered."); for (std::list::iterator it = synchronizers.begin(); it != synchronizers.end(); ++it) { (*it)->asynchronousSynchronize(tag); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void GhostSynchronizer::waitEndSynchronize(GhostSynchronizationTag tag) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(registered_synchronization.find(tag) != registered_synchronization.end(), "Tag " << tag << " is not registered."); for (std::list::iterator it = synchronizers.begin(); it != synchronizers.end(); ++it) { (*it)->waitEndSynchronize(tag); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void GhostSynchronizer::registerTag(GhostSynchronizationTag tag, const std::string & name) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(registered_synchronization.find(tag) == registered_synchronization.end(), "Tag " << tag << " (" << name << ") already registered."); registered_synchronization[tag] = name; for (std::list::iterator it = synchronizers.begin(); it != synchronizers.end(); ++it) { (*it)->registerTag(tag); } AKANTU_DEBUG_OUT(); } +/* -------------------------------------------------------------------------- */ +void GhostSynchronizer::allReduce(Real * values, + const SynchronizerOperation & op, + UInt nb_values) { + AKANTU_DEBUG_IN(); + + for (std::list::iterator it = synchronizers.begin(); + it != synchronizers.end(); + ++it) { + (*it)->allReduce(values, nb_values, op); + } + + AKANTU_DEBUG_OUT(); +} + /* -------------------------------------------------------------------------- */ void GhostSynchronizer::registerSynchronizer(Synchronizer & synchronizer) { AKANTU_DEBUG_IN(); synchronizers.push_back(&synchronizer); synchronizer.registerGhostSynchronizer(*this); std::map::iterator it; for (it = registered_synchronization.begin(); it != registered_synchronization.end(); ++it) { synchronizer.registerTag(it->first); } AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/synchronizer/ghost_synchronizer.hh b/synchronizer/ghost_synchronizer.hh index bf91dd1bb..cceb2e94b 100644 --- a/synchronizer/ghost_synchronizer.hh +++ b/synchronizer/ghost_synchronizer.hh @@ -1,111 +1,116 @@ /** * @file ghost_synchronizer.hh * @author Nicolas Richart * @date Fri Aug 20 17:40:08 2010 * * @brief Class of ghost synchronisation (PBC or parallel communication) * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_GHOST_SYNCHRONIZER_HH__ #define __AKANTU_GHOST_SYNCHRONIZER_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ namespace akantu { class Synchronizer; } __BEGIN_AKANTU__ class GhostSynchronizer { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: GhostSynchronizer(); virtual ~GhostSynchronizer(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: + /* ------------------------------------------------------------------------ */ /// synchronize ghosts void synchronize(GhostSynchronizationTag tag); /// asynchronous synchronization of ghosts void asynchronousSynchronize(GhostSynchronizationTag tag); /// wait end of asynchronous synchronization of ghosts void waitEndSynchronize(GhostSynchronizationTag tag); + /// reduce a value (essentially for communications synchronizer) + void allReduce(Real * values, const SynchronizerOperation & op, UInt nb_values = 1); + + /* ------------------------------------------------------------------------ */ /// register a new communication void registerTag(GhostSynchronizationTag tag, const std::string & name); /// register a new synchronization void registerSynchronizer(Synchronizer & synchronizer); public: virtual UInt getNbDataToPack(const Element & element, GhostSynchronizationTag tag) const = 0; virtual UInt getNbDataToUnpack(const Element & element, GhostSynchronizationTag tag) const = 0; virtual void packData(Real ** buffer, const Element & element, GhostSynchronizationTag tag) const = 0; virtual void unpackData(Real ** buffer, const Element & element, GhostSynchronizationTag tag) const = 0; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// list of synchronizers std::list synchronizers; /// number of tags registered UInt nb_ghost_synchronization_tags; /// list of registered synchronization std::map registered_synchronization; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "ghost_synchronizer_inline_impl.cc" /// standard output stream operator // inline std::ostream & operator <<(std::ostream & stream, const GhostSynchronizer & _this) // { // _this.printself(stream); // return stream; // } __END_AKANTU__ #endif /* __AKANTU_GHOST_SYNCHRONIZER_HH__ */ diff --git a/synchronizer/static_communicator.hh b/synchronizer/static_communicator.hh index 4e9601e8e..6a07bbf39 100644 --- a/synchronizer/static_communicator.hh +++ b/synchronizer/static_communicator.hh @@ -1,129 +1,132 @@ /** * @file static_communicator.hh * @author Nicolas Richart * @date Thu Aug 19 15:34:09 2010 * * @brief Class handling the parallel communications * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_STATIC_COMMUNICATOR_HH__ #define __AKANTU_STATIC_COMMUNICATOR_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ class CommunicationRequest { public: CommunicationRequest(UInt source, UInt dest); virtual ~CommunicationRequest(); virtual void printself(std::ostream & stream, int indent = 0) const; AKANTU_GET_MACRO(Source, source, UInt); AKANTU_GET_MACRO(Destination, destination, UInt); private: UInt source; UInt destination; UInt id; static UInt counter; }; class StaticCommunicator { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ protected: StaticCommunicator() { }; public: virtual ~StaticCommunicator() { }; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual void send(UInt * buffer, Int size, Int receiver, Int tag) = 0; virtual void send(Real * buffer, Int size, Int receiver, Int tag) = 0; virtual void receive(UInt * buffer, Int size, Int sender, Int tag) = 0; virtual void receive(Real * buffer, Int size, Int sender, Int tag) = 0; virtual CommunicationRequest * asyncSend(UInt * buffer, Int size, Int receiver, Int tag) = 0; virtual CommunicationRequest * asyncSend(Real * buffer, Int size, Int receiver, Int tag) = 0; virtual CommunicationRequest * asyncReceive(UInt * buffer, Int size, Int sender, Int tag) = 0; virtual CommunicationRequest * asyncReceive(Real * buffer, Int size, Int sender, Int tag) = 0; virtual bool testRequest(CommunicationRequest * request) = 0; virtual void wait(CommunicationRequest * request) = 0; virtual void waitAll(std::vector & requests) = 0; virtual void freeCommunicationRequest(CommunicationRequest * request); virtual void freeCommunicationRequest(std::vector & requests); virtual void barrier() = 0; + virtual void allReduce(Real * values, UInt nb_values, const SynchronizerOperation & op) = 0; + virtual void allReduce(UInt * values, UInt nb_values, const SynchronizerOperation & op) = 0; + /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: virtual Int getNbProc() = 0; virtual Int whoAmI() = 0; static StaticCommunicator * getStaticCommunicator(CommunicatorType type = _communicator_mpi); static StaticCommunicator * getStaticCommunicator(int * argc, char *** argv, CommunicatorType type = _communicator_mpi); static bool isInstantiated() { return is_instantiated; }; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: static bool is_instantiated; static StaticCommunicator * static_communicator; protected: Int prank; Int psize; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "static_communicator_inline_impl.cc" /* -------------------------------------------------------------------------- */ /* Inline Functions VectorBase */ /* -------------------------------------------------------------------------- */ inline std::ostream & operator<<(std::ostream & stream, const CommunicationRequest & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_STATIC_COMMUNICATOR_HH__ */ diff --git a/synchronizer/static_communicator_dummy.hh b/synchronizer/static_communicator_dummy.hh index bd2a0346b..816891e31 100644 --- a/synchronizer/static_communicator_dummy.hh +++ b/synchronizer/static_communicator_dummy.hh @@ -1,86 +1,89 @@ /** * @file static_communicator_dummy.hh * @author Nicolas Richart * @date Thu Aug 19 15:34:09 2010 * * @brief Class handling the parallel communications * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_STATIC_COMMUNICATOR_DUMMY_HH__ #define __AKANTU_STATIC_COMMUNICATOR_DUMMY_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ class StaticCommunicatorDummy : public StaticCommunicator { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: virtual ~StaticCommunicatorDummy() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual void send(UInt * buffer, Int size, Int receiver, Int tag) {}; virtual void send(Real * buffer, Int size, Int receiver, Int tag) {}; virtual void receive(UInt * buffer, Int size, Int sender, Int tag) {}; virtual void receive(Real * buffer, Int size, Int sender, Int tag) {}; virtual CommunicationRequest * asyncSend(UInt * buffer, Int size, Int receiver, Int tag) { return new CommunicationRequest(0, 0); }; virtual CommunicationRequest * asyncSend(Real * buffer, Int size, Int receiver, Int tag) { return new CommunicationRequest(0, 0); }; virtual CommunicationRequest * asyncReceive(UInt * buffer, Int size, Int sender, Int tag) { return new CommunicationRequest(0, 0); }; virtual CommunicationRequest * asyncReceive(Real * buffer, Int size, Int sender, Int tag) { return new CommunicationRequest(0, 0); }; virtual bool testRequest(CommunicationRequest * request) { return true; }; virtual void wait(CommunicationRequest * request) {}; virtual void waitAll(std::vector & requests) {}; virtual void barrier() {}; + virtual void allReduce(Real * values, UInt nb_values, const SynchronizerOperation & op) {}; + virtual void allReduce(UInt * values, UInt nb_values, const SynchronizerOperation & op) {}; + /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: virtual Int getNbProc() { return 1; }; virtual Int whoAmI() { return 0; }; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ }; __END_AKANTU__ #endif /* __AKANTU_STATIC_COMMUNICATOR_DUMMY_HH__ */ diff --git a/synchronizer/static_communicator_mpi.cc b/synchronizer/static_communicator_mpi.cc new file mode 100644 index 000000000..435178454 --- /dev/null +++ b/synchronizer/static_communicator_mpi.cc @@ -0,0 +1,28 @@ +/** + * @file static_communicator_mpi.cc + * @author Nicolas Richart + * @date Tue Sep 14 11:37:10 2010 + * + * @brief + * + * @section LICENSE + * + * + * + */ + +/* -------------------------------------------------------------------------- */ +#include "static_communicator_mpi.hh" + +/* -------------------------------------------------------------------------- */ + +__BEGIN_AKANTU__ + +MPI_Op StaticCommunicatorMPI::synchronizer_operation_to_mpi_op[_so_null + 1] = { + MPI_SUM, + MPI_MIN, + MPI_MAX, + MPI_OP_NULL +}; + +__END_AKANTU__ diff --git a/synchronizer/static_communicator_mpi.hh b/synchronizer/static_communicator_mpi.hh index e39a781c9..2db610ada 100644 --- a/synchronizer/static_communicator_mpi.hh +++ b/synchronizer/static_communicator_mpi.hh @@ -1,105 +1,110 @@ /** * @file static_communicator_mpi.hh * @author Nicolas Richart * @date Thu Sep 2 19:59:58 2010 * * @brief class handling parallel communication trough MPI * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_STATIC_COMMUNICATOR_MPI_HH__ #define __AKANTU_STATIC_COMMUNICATOR_MPI_HH__ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include "static_communicator.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ class CommunicationRequestMPI : public CommunicationRequest { public: inline CommunicationRequestMPI(UInt source, UInt dest); inline ~CommunicationRequestMPI(); inline MPI_Request * getMPIRequest() { return request; }; private: MPI_Request * request; }; /* -------------------------------------------------------------------------- */ class StaticCommunicatorMPI : public virtual StaticCommunicator { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: inline StaticCommunicatorMPI(int * argc, char *** argv); inline virtual ~StaticCommunicatorMPI(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: inline void send(UInt * buffer, Int size, Int receiver, Int tag); inline void send(Real * buffer, Int size, Int receiver, Int tag); inline void receive(UInt * buffer, Int size, Int sender, Int tag); inline void receive(Real * buffer, Int size, Int sender, Int tag); inline CommunicationRequest * asyncSend(UInt * buffer, Int size, Int receiver, Int tag); inline CommunicationRequest * asyncSend(Real * buffer, Int size, Int receiver, Int tag); - virtual CommunicationRequest * asyncReceive(UInt * buffer, Int size, + inline CommunicationRequest * asyncReceive(UInt * buffer, Int size, Int sender, Int tag); - virtual CommunicationRequest * asyncReceive(Real * buffer, Int size, + inline CommunicationRequest * asyncReceive(Real * buffer, Int size, Int sender, Int tag); - virtual bool testRequest(CommunicationRequest * request); + inline bool testRequest(CommunicationRequest * request); inline void wait(CommunicationRequest * request); inline void waitAll(std::vector & requests); inline void barrier(); + inline void allReduce(Real * values, UInt nb_values, const SynchronizerOperation & op); + inline void allReduce(UInt * values, UInt nb_values, const SynchronizerOperation & op); + /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: inline void setMPIComm(MPI_Comm comm); inline Int getNbProc() { return psize; }; inline Int whoAmI() { return prank; }; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: MPI_Comm communicator; + + static MPI_Op synchronizer_operation_to_mpi_op[_so_null + 1]; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "static_communicator_mpi_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_STATIC_COMMUNICATOR_MPI_HH__ */ diff --git a/synchronizer/static_communicator_mpi_inline_impl.cc b/synchronizer/static_communicator_mpi_inline_impl.cc index df3b558bc..78e862342 100644 --- a/synchronizer/static_communicator_mpi_inline_impl.cc +++ b/synchronizer/static_communicator_mpi_inline_impl.cc @@ -1,143 +1,163 @@ /** * @file static_communicator_mpi_inline_impl.hh * @author Nicolas Richart * @date Thu Sep 2 15:10:51 2010 * * @brief implementation of the mpi communicator * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ inline CommunicationRequestMPI::CommunicationRequestMPI(UInt source, UInt dest) : CommunicationRequest(source, dest) { request = new MPI_Request; } /* -------------------------------------------------------------------------- */ inline CommunicationRequestMPI::~CommunicationRequestMPI() { delete request; } /* -------------------------------------------------------------------------- */ inline StaticCommunicatorMPI::StaticCommunicatorMPI(int * argc, char *** argv) { MPI_Init(argc, argv); setMPIComm(MPI_COMM_WORLD); }; /* -------------------------------------------------------------------------- */ inline StaticCommunicatorMPI::~StaticCommunicatorMPI() { MPI_Finalize(); }; /* -------------------------------------------------------------------------- */ inline void StaticCommunicatorMPI::setMPIComm(MPI_Comm comm) { communicator = comm; MPI_Comm_rank(communicator, &prank); MPI_Comm_size(communicator, &psize); } /* -------------------------------------------------------------------------- */ inline void StaticCommunicatorMPI::send(UInt * buffer, Int size, Int receiver, Int tag) { - Int ret = MPI_Send(buffer, size, MPI_UNSIGNED, receiver, tag, communicator); + int ret = MPI_Send(buffer, size, MPI_UNSIGNED, receiver, tag, communicator); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Send."); } /* -------------------------------------------------------------------------- */ inline void StaticCommunicatorMPI::send(Real * buffer, Int size, Int receiver, Int tag) { - Int ret = MPI_Send(buffer, size, MPI_DOUBLE, receiver, tag, communicator); + int ret = MPI_Send(buffer, size, MPI_DOUBLE, receiver, tag, communicator); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Send."); } /* -------------------------------------------------------------------------- */ inline void StaticCommunicatorMPI::receive(UInt * buffer, Int size, Int sender, Int tag) { MPI_Status status; - Int ret = MPI_Recv(buffer, size, MPI_UNSIGNED, sender, tag, communicator, &status); + int ret = MPI_Recv(buffer, size, MPI_UNSIGNED, sender, tag, communicator, &status); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Recv."); } /* -------------------------------------------------------------------------- */ inline void StaticCommunicatorMPI::receive(Real * buffer, Int size, Int sender, Int tag) { MPI_Status status; - Int ret = MPI_Recv(buffer, size, MPI_DOUBLE, sender, tag, communicator, &status); + int ret = MPI_Recv(buffer, size, MPI_DOUBLE, sender, tag, communicator, &status); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Recv."); } /* -------------------------------------------------------------------------- */ inline CommunicationRequest * StaticCommunicatorMPI::asyncSend(UInt * buffer, Int size, Int receiver, Int tag) { CommunicationRequestMPI * request = new CommunicationRequestMPI(prank, receiver); - Int ret = MPI_Isend(buffer, size, MPI_UNSIGNED, receiver, tag, communicator, request->getMPIRequest()); + int ret = MPI_Isend(buffer, size, MPI_UNSIGNED, receiver, tag, communicator, request->getMPIRequest()); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Isend."); return request; }; /* -------------------------------------------------------------------------- */ inline CommunicationRequest * StaticCommunicatorMPI::asyncSend(Real * buffer, Int size, Int receiver, Int tag) { CommunicationRequestMPI * request = new CommunicationRequestMPI(prank, receiver); - Int ret = MPI_Isend(buffer, size, MPI_DOUBLE, receiver, tag, communicator, request->getMPIRequest()); + int ret = MPI_Isend(buffer, size, MPI_DOUBLE, receiver, tag, communicator, request->getMPIRequest()); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Isend."); return request; }; /* -------------------------------------------------------------------------- */ inline CommunicationRequest * StaticCommunicatorMPI::asyncReceive(UInt * buffer, Int size, Int sender, Int tag) { CommunicationRequestMPI * request = new CommunicationRequestMPI(sender, prank); - Int ret = MPI_Irecv(buffer, size, MPI_UNSIGNED, sender, tag, communicator, request->getMPIRequest()); + int ret = MPI_Irecv(buffer, size, MPI_UNSIGNED, sender, tag, communicator, request->getMPIRequest()); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Irecv."); return request; }; /* -------------------------------------------------------------------------- */ inline CommunicationRequest * StaticCommunicatorMPI::asyncReceive(Real * buffer, Int size, Int sender, Int tag) { CommunicationRequestMPI * request = new CommunicationRequestMPI(sender, prank); - Int ret = MPI_Irecv(buffer, size, MPI_DOUBLE, sender, tag, communicator, request->getMPIRequest()); + int ret = MPI_Irecv(buffer, size, MPI_DOUBLE, sender, tag, communicator, request->getMPIRequest()); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Irecv."); return request; }; /* -------------------------------------------------------------------------- */ inline bool StaticCommunicatorMPI::testRequest(CommunicationRequest * request) { MPI_Status status; + int flag; CommunicationRequestMPI * req_mpi = static_cast(request); MPI_Request * req = req_mpi->getMPIRequest(); - int flag; - Int ret = MPI_Test(req, &flag, &status); + + int ret = MPI_Test(req, &flag, &status); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Test."); return (flag != 0); } /* -------------------------------------------------------------------------- */ inline void StaticCommunicatorMPI::wait(CommunicationRequest * request) { MPI_Status status; CommunicationRequestMPI * req_mpi = static_cast(request); MPI_Request * req = req_mpi->getMPIRequest(); - Int ret = MPI_Wait(req, &status); + + int ret = MPI_Wait(req, &status); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Wait."); } /* -------------------------------------------------------------------------- */ inline void StaticCommunicatorMPI::waitAll(std::vector & requests) { MPI_Status status; std::vector::iterator it; for(it = requests.begin(); it != requests.end(); ++it) { MPI_Request * req = static_cast(*it)->getMPIRequest(); - Int ret = MPI_Wait(req, &status); + int ret = MPI_Wait(req, &status); AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Wait."); } } /* -------------------------------------------------------------------------- */ inline void StaticCommunicatorMPI::barrier() { MPI_Barrier(communicator); } + +/* -------------------------------------------------------------------------- */ +inline void StaticCommunicatorMPI::allReduce(UInt * values, UInt nb_values, + const SynchronizerOperation & op) { + int ret = MPI_Allreduce(MPI_IN_PLACE, values, nb_values, MPI_UNSIGNED, + synchronizer_operation_to_mpi_op[op], + communicator); + AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Allreduce."); +} + +/* -------------------------------------------------------------------------- */ +inline void StaticCommunicatorMPI::allReduce(Real * values, UInt nb_values, + const SynchronizerOperation & op) { + int ret = MPI_Allreduce(MPI_IN_PLACE, values, nb_values, MPI_DOUBLE, + synchronizer_operation_to_mpi_op[op], + communicator); + AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Allreduce."); +} diff --git a/synchronizer/synchronizer.hh b/synchronizer/synchronizer.hh index 0d08a4ac7..c2f8f0a1f 100644 --- a/synchronizer/synchronizer.hh +++ b/synchronizer/synchronizer.hh @@ -1,79 +1,82 @@ /** * @file synchronizer.hh * @author Nicolas Richart * @date Mon Aug 23 13:48:37 2010 * * @brief interface for communicator and pbc synchronizers * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SYNCHRONIZER_HH__ #define __AKANTU_SYNCHRONIZER_HH__ /* -------------------------------------------------------------------------- */ #include "aka_memory.hh" /* -------------------------------------------------------------------------- */ namespace akantu { class GhostSynchronizer; }; __BEGIN_AKANTU__ class Synchronizer : public Memory { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Synchronizer(SynchronizerID id = "synchronizer", MemoryID memory_id = 0); virtual ~Synchronizer() { }; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// synchronize ghosts virtual void synchronize(GhostSynchronizationTag tag) = 0; /// asynchronous synchronization of ghosts virtual void asynchronousSynchronize(GhostSynchronizationTag tag) = 0; /// wait end of asynchronous synchronization of ghosts virtual void waitEndSynchronize(GhostSynchronizationTag tag) = 0; + virtual void allReduce(Real * values, UInt nb_values, const SynchronizerOperation & op) = 0; + + /* ------------------------------------------------------------------------ */ /// register a new communication virtual void registerTag(GhostSynchronizationTag tag) = 0; void registerGhostSynchronizer(const GhostSynchronizer & ghost_synchronizer); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id of the synchronizer SynchronizerID id; /// the associated ghost_synchonizer const GhostSynchronizer * ghost_synchronizer; }; __END_AKANTU__ #endif /* __AKANTU_SYNCHRONIZER_HH__ */ diff --git a/test/test_fem/CMakeLists.txt b/test/test_fem/CMakeLists.txt index 1c1a7e034..0b8c4ebe7 100644 --- a/test/test_fem/CMakeLists.txt +++ b/test/test_fem/CMakeLists.txt @@ -1,57 +1,57 @@ #=============================================================================== # @file CMakeLists.txt # @author Guillaume Anciaux # @date Fri Jun 11 09:46:59 2010 # # @section LICENSE # # # # @section DESCRIPTION # #=============================================================================== #=============================================================================== # 1D #=============================================================================== add_mesh(test_fem_line1_mesh line.geo 1 1) -REGISTER_TEST(test_integrate_line_1 test_integrate_line_1.cc) +register_test(test_integrate_line_1 test_integrate_line_1.cc) add_dependencies(test_integrate_line_1 line) #=============================================================================== register_test(test_interpolate_line_1 test_interpolate_line_1.cc) add_dependencies(test_interpolate_line_1 test_fem_line1_mesh) #=============================================================================== register_test(test_gradient_line_1 test_gradient_line_1.cc) add_dependencies(test_gradient_line_1 test_fem_line1_mesh) -REGISTER_TEST(test_integrate_triangle_1 test_integrate_triangle_1.cc) +register_test(test_integrate_triangle_1 test_integrate_triangle_1.cc) add_dependencies(test_integrate_triangle_1 ltriangle) #=============================================================================== # 2D #=============================================================================== add_mesh(test_fem_triangle1_mesh triangle.geo 2 1) -REGISTER_TEST(test_integrate_tetrahedra_1 test_integrate_tetrahedra_1.cc) +register_test(test_integrate_tetrahedra_1 test_integrate_tetrahedra_1.cc) add_dependencies(test_integrate_tetrahedra_1 llcube) #=============================================================================== register_test(test_interpolate_triangle_1 test_interpolate_triangle_1.cc) add_dependencies(test_interpolate_triangle_1 test_fem_triangle1_mesh) #=============================================================================== register_test(test_gradient_triangle_1 test_gradient_triangle_1.cc) add_dependencies(test_gradient_triangle_1 test_fem_triangle1_mesh) #=============================================================================== # 3D #=============================================================================== add_mesh(test_fem_tetra1_mesh cube.geo 3 1) #=============================================================================== register_test(test_interpolate_tetrahedra_1 test_interpolate_tetrahedra_1.cc) add_dependencies(test_interpolate_tetrahedra_1 test_fem_tetra1_mesh) #=============================================================================== register_test(test_gradient_tetrahedra_1 test_gradient_tetrahedra_1.cc) add_dependencies(test_gradient_tetrahedra_1 test_fem_tetra1_mesh) diff --git a/test/test_solid_mechanics_model/CMakeLists.txt b/test/test_solid_mechanics_model/CMakeLists.txt index 3bdf40e74..bc5012121 100644 --- a/test/test_solid_mechanics_model/CMakeLists.txt +++ b/test/test_solid_mechanics_model/CMakeLists.txt @@ -1,50 +1,58 @@ #=============================================================================== # @file CMakeLists.txt # @author Guillaume Anciaux # @date Fri Jun 11 09:46:59 2010 # # @section LICENSE # # # # @section DESCRIPTION # #=============================================================================== #=============================================================================== add_mesh(test_solid_mechanics_model_square_mesh triangle.geo 2 1) add_mesh(test_solid_mechanics_model_circle_mesh circle.geo 2 1) register_test(test_solid_mechanics_model_square test_solid_mechanics_model_square.cc) add_dependencies(test_solid_mechanics_model_square test_solid_mechanics_model_square_mesh test_solid_mechanics_model_circle_mesh) #=============================================================================== add_mesh(test_bar_traction_2d_mesh bar.geo 2 1) register_test(test_solid_mechanics_model_bar_traction2d test_solid_mechanics_model_bar_traction2d.cc) add_dependencies(test_solid_mechanics_model_bar_traction2d test_bar_traction_2d_mesh) if(AKANTU_SCOTCH_ON AND AKANTU_MPI_ON) register_test(test_solid_mechanics_model_bar_traction2d_parallel test_solid_mechanics_model_bar_traction2d_parallel.cc) add_dependencies(test_solid_mechanics_model_bar_traction2d_parallel test_bar_traction_2d_mesh) + + + add_mesh(test_solid_mechanics_model_line_mesh line.geo 1 1) + + register_test(test_solid_mechanics_model_line_parallel + test_solid_mechanics_model_line_parallel.cc) + add_dependencies(test_solid_mechanics_model_line_parallel + test_solid_mechanics_model_line_mesh) endif() #=============================================================================== add_mesh(test_cube3d_mesh cube.geo 3 1) register_test(test_solid_mechanics_model_cube3d test_solid_mechanics_model_cube3d.cc) add_dependencies(test_solid_mechanics_model_cube3d test_cube3d_mesh) #=============================================================================== configure_file( ${CMAKE_CURRENT_SOURCE_DIR}/material.dat ${CMAKE_CURRENT_BINARY_DIR}/material.dat COPYONLY ) \ No newline at end of file diff --git a/test/test_solid_mechanics_model/bar.geo b/test/test_solid_mechanics_model/bar.geo index 1c0161ae3..fa76fdd0c 100644 --- a/test/test_solid_mechanics_model/bar.geo +++ b/test/test_solid_mechanics_model/bar.geo @@ -1,42 +1,30 @@ -// CUBE ON CUBE (2 BODIES) 3D: - -// REMARKS: -// Physical Surfaces are defined so that load_mesh_msh can load -// directly the surface element structure. When creating Plane Surface, -// the surface normal has to point to the inside of the body. Otherwise the -// solvecontact3d algorithm won't work. The Physical Surface number -// corresponds to the face number. The face numbering has to start with 1! -// Physical Volume defines the material of each body. This can be read -// by the load_mesh_msh function in adlib. - // Mesh size -h = 0.05; // Top cube -hy = 0.5; // Top cube +h = 0.1; // Top cube // Dimensions of top cube Lx = 10; Ly = 1; // ------------------------------------------ // Geometry // ------------------------------------------ // Base Cube Point(101) = { 0.0, 0.0, 0.0, h}; // Bottom Face Point(102) = { Lx, 0.0, 0.0, h}; // Bottom Face Point(103) = { Lx, Ly, 0.0, h}; // Bottom Face Point(104) = { 0.0, Ly, 0.0, h}; // Bottom Face // Base Cube Line(101) = {101,102}; // Bottom Face Line(102) = {102,103}; // Bottom Face Line(103) = {103,104}; // Bottom Face Line(104) = {104,101}; // Bottom Face // Base Cube Line Loop(101) = {101:104}; Plane Surface(101) = {101}; Physical Surface(1) = {101}; diff --git a/test/test_solid_mechanics_model/line.geo b/test/test_solid_mechanics_model/line.geo new file mode 100644 index 000000000..056ecc24a --- /dev/null +++ b/test/test_solid_mechanics_model/line.geo @@ -0,0 +1,9 @@ +// Mesh size +h = 0.01; + +Lx = 1; + +Point(1) = { 0.0, 0.0, 0.0, h}; +Point(2) = { Lx, 0.0, 0.0, h}; + +Line(1) = {1,2}; diff --git a/test/test_solid_mechanics_model/test_solid_mechanics_model_bar_traction2d_parallel.cc b/test/test_solid_mechanics_model/test_solid_mechanics_model_bar_traction2d_parallel.cc index f6d133faf..9e470e5d9 100644 --- a/test/test_solid_mechanics_model/test_solid_mechanics_model_bar_traction2d_parallel.cc +++ b/test/test_solid_mechanics_model/test_solid_mechanics_model_bar_traction2d_parallel.cc @@ -1,228 +1,218 @@ /** * @file test_solid_mechanics_model.cc * @author Nicolas Richart * @date Tue Jul 27 14:34:13 2010 * * @brief test of the class SolidMechanicsModel * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "static_communicator.hh" #include "communicator.hh" #include "mesh_partition_scotch.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER #define CHECK_STRESS - int main(int argc, char *argv[]) { + akantu::ElementType type = akantu::_triangle_1; akantu::UInt spatial_dimension = 2; akantu::UInt max_steps = 1000; akantu::Real time_factor = 0.2; akantu::initialize(&argc, &argv); - akantu::debug::setDebugLevel(akantu::dblDump); + akantu::debug::setDebugLevel(akantu::dblWarning); // akantu::Real epot, ekin; akantu::Mesh mesh(spatial_dimension); akantu::StaticCommunicator * comm = akantu::StaticCommunicator::getStaticCommunicator(); akantu::Int psize = comm->getNbProc(); akantu::Int prank = comm->whoAmI(); akantu::Communicator * communicator; if(prank == 0) { akantu::MeshIOMSH mesh_io; mesh_io.read("bar.msh", mesh); - akantu::MeshPartition * partition = + akantu::MeshPartition * partition = new akantu::MeshPartitionScotch(mesh, spatial_dimension); partition->partitionate(psize); communicator = akantu::Communicator::createCommunicatorDistributeMesh(mesh, partition); delete partition; } else { communicator = akantu::Communicator::createCommunicatorDistributeMesh(mesh, NULL); } + // std::cout << mesh << std::endl; + akantu::SolidMechanicsModel * model = new akantu::SolidMechanicsModel(mesh); akantu::UInt nb_nodes = model->getFEM().getMesh().getNbNodes(); - akantu::UInt nb_element = model->getFEM().getMesh().getNbElement(akantu::_triangle_1); + akantu::UInt nb_element = model->getFEM().getMesh().getNbElement(type); /// model initialization model->initVectors(); /// set vectors to 0 memset(model->getForce().values, 0, spatial_dimension*nb_nodes*sizeof(akantu::Real)); memset(model->getVelocity().values, 0, spatial_dimension*nb_nodes*sizeof(akantu::Real)); memset(model->getAcceleration().values, 0, spatial_dimension*nb_nodes*sizeof(akantu::Real)); memset(model->getDisplacement().values, 0, spatial_dimension*nb_nodes*sizeof(akantu::Real)); model->readMaterials("material.dat"); model->initMaterials(); model->registerSynchronizer(*communicator); model->initModel(); - std::cout << model->getMaterial(0) << std::endl; + if(prank == 0) + std::cout << model->getMaterial(0) << std::endl; model->assembleMass(); -#ifdef AKANTU_USE_IOHELPER - /// set to 0 only for the first paraview dump - memset(model->getResidual().values, 0, - spatial_dimension*nb_nodes*sizeof(akantu::Real)); - memset(model->getMaterial(0).getStrain(akantu::_triangle_1).values, 0, - spatial_dimension*spatial_dimension*nb_element*sizeof(akantu::Real)); - memset(model->getMaterial(0).getStress(akantu::_triangle_1).values, 0, - spatial_dimension*spatial_dimension*nb_element*sizeof(akantu::Real)); -#endif //AKANTU_USE_IOHELPER - - /// boundary conditions akantu::Real eps = 1e-16; for (akantu::UInt i = 0; i < nb_nodes; ++i) { if(model->getFEM().getMesh().getNodes().values[spatial_dimension*i] >= 9) - model->getDisplacement().values[spatial_dimension*i] = (model->getFEM().getMesh().getNodes().values[spatial_dimension*i] - 9) / 100. ; + model->getDisplacement().values[spatial_dimension*i] + = (model->getFEM().getMesh().getNodes().values[spatial_dimension*i] - 9) / 100.; if(model->getFEM().getMesh().getNodes().values[spatial_dimension*i] <= eps) - model->getBoundary().values[spatial_dimension*i] = true; + model->getBoundary().values[spatial_dimension*i] = true; if(model->getFEM().getMesh().getNodes().values[spatial_dimension*i + 1] <= eps || model->getFEM().getMesh().getNodes().values[spatial_dimension*i + 1] >= 1 - eps ) { model->getBoundary().values[spatial_dimension*i + 1] = true; } } + model->synchronizeBoundaries(); + akantu::Real time_step = model->getStableTimeStep() * time_factor; - std::cout << "Time Step = " << time_step << "s" << std::endl; + if(prank == 0) + std::cout << "Time Step = " << time_step << "s" << std::endl; model->setTimeStep(time_step); - // model->setTimeStep(3.54379e-07); #ifdef AKANTU_USE_IOHELPER DumperParaview dumper; dumper.SetMode(TEXT); + dumper.SetParallelContext(prank, psize); dumper.SetPoints(model->getFEM().getMesh().getNodes().values, - spatial_dimension, nb_nodes, "coordinates_par"); - dumper.SetConnectivity((int *)model->getFEM().getMesh().getConnectivity(akantu::_triangle_1).values, + spatial_dimension, nb_nodes, "bar2d_para"); + dumper.SetConnectivity((int *)model->getFEM().getMesh().getConnectivity(type).values, TRIANGLE1, nb_element, C_MODE); dumper.AddNodeDataField(model->getDisplacement().values, spatial_dimension, "displacements"); + dumper.AddNodeDataField(model->getMass().values, + 1, "mass"); dumper.AddNodeDataField(model->getVelocity().values, spatial_dimension, "velocity"); dumper.AddNodeDataField(model->getResidual().values, spatial_dimension, "force"); - dumper.AddElemDataField(model->getMaterial(0).getStrain(akantu::_triangle_1).values, + dumper.AddNodeDataField(model->getAcceleration().values, + spatial_dimension, "acceleration"); + dumper.AddElemDataField(model->getMaterial(0).getStrain(type).values, spatial_dimension*spatial_dimension, "strain"); - dumper.AddElemDataField(model->getMaterial(0).getStress(akantu::_triangle_1).values, + dumper.AddElemDataField(model->getMaterial(0).getStress(type).values, spatial_dimension*spatial_dimension, "stress"); + double * part = new double[nb_element]; + for (unsigned int i = 0; i < nb_element; ++i) + part[i] = prank; + dumper.AddElemDataField(part, 1, "partitions"); dumper.SetEmbeddedValue("displacements", 1); dumper.SetPrefix("paraview/"); dumper.Init(); - dumper.Dump(); + //dumper.Dump(); + + DumperParaview dumper_ghost; + if (psize > 1) { + unsigned int nb_ghost_element = mesh.getNbGhostElement(type); + dumper_ghost.SetMode(TEXT); + dumper_ghost.SetParallelContext(prank, psize); + dumper_ghost.SetPoints(mesh.getNodes().values, + spatial_dimension, nb_nodes, "bar2d_para_ghost"); + dumper_ghost.SetConnectivity((int*) mesh.getGhostConnectivity(type).values, + TRIANGLE1, nb_ghost_element, C_MODE); + dumper_ghost.AddNodeDataField(model->getDisplacement().values, + spatial_dimension, "displacements"); + dumper_ghost.AddNodeDataField(model->getMass().values, + 1, "mass"); + dumper_ghost.AddNodeDataField(model->getVelocity().values, + spatial_dimension, "velocity"); + dumper_ghost.AddNodeDataField(model->getResidual().values, + spatial_dimension, "force"); + dumper_ghost.AddNodeDataField(model->getAcceleration().values, + spatial_dimension, "acceleration"); + dumper_ghost.SetPrefix("paraview/"); + dumper_ghost.Init(); + // dumper_ghost.Dump(); + } #endif //AKANTU_USE_IOHELPER + model->setPotentialEnergyFlagOn(); -#ifdef CHECK_STRESS - std::ofstream outfile; - outfile.open("stress"); -#endif // CHECK_STRESS + std::ofstream energy; + if(prank == 0) { + energy.open("energy.csv"); + energy << "id,epot,ekin,tot" << std::endl; + } - model->setPotentialEnergyFlagOn(); for(akantu::UInt s = 1; s <= max_steps; ++s) { model->explicitPred(); model->updateResidual(); model->updateAcceleration(); model->explicitCorr(); - // epot = model->getPotentialEnergy(); - // ekin = model->getKineticEnergy(); - // std::cout << s << " " << epot << " " << ekin << " " << epot + ekin - // << std::endl; - -#ifdef CHECK_STRESS - akantu::Real max_stress = std::numeric_limits::min(); - akantu::UInt max_el = 0; - akantu::Real * stress = model->getMaterial(0).getStress(akantu::_triangle_1).values; - for (akantu::UInt i = 0; i < nb_element; ++i) { - if(max_stress < stress[i*spatial_dimension*spatial_dimension]) { - max_stress = stress[i*spatial_dimension*spatial_dimension]; - max_el = i; - } - } + akantu::Real epot = model->getPotentialEnergy(); + akantu::Real ekin = model->getKineticEnergy(); - akantu::Real * coord = model->getFEM().getMesh().getNodes().values; - akantu::Real * disp_val = model->getDisplacement().values; - akantu::UInt * conn = model->getFEM().getMesh().getConnectivity(akantu::_triangle_1).values; - akantu::UInt nb_nodes_per_element = model->getFEM().getMesh().getNbNodesPerElement(akantu::_triangle_1); - akantu::Real * coords = new akantu::Real[spatial_dimension]; - akantu::Real min_x = std::numeric_limits::max(); - akantu::Real max_x = std::numeric_limits::min(); - - akantu::Real stress_range = 5e7; - for (akantu::UInt el = 0; el < nb_element; ++el) { - if(stress[el*spatial_dimension*spatial_dimension] > max_stress - stress_range) { - akantu::UInt el_offset = el * nb_nodes_per_element; - memset(coords, 0, spatial_dimension*sizeof(akantu::Real)); - for (akantu::UInt n = 0; n < nb_nodes_per_element; ++n) { - for (akantu::UInt i = 0; i < spatial_dimension; ++i) { - akantu::UInt node = conn[el_offset + n] * spatial_dimension; - coords[i] += (coord[node + i] + disp_val[node + i]) - / ((akantu::Real) nb_nodes_per_element); - } - } - min_x = min_x < coords[0] ? min_x : coords[0]; - max_x = max_x > coords[0] ? max_x : coords[0]; - } + if(prank == 0 && (s % 10 == 0)) { + energy << s << "," << epot << "," << ekin << "," << epot + ekin + << std::endl; } - - - outfile << s << " " << .5 * (min_x + max_x) << " " << min_x << " " << max_x << " " << max_x - min_x << " " << max_stress << std::endl; - - delete [] coords; -#endif // CHECK_STRESS - - #ifdef AKANTU_USE_IOHELPER - if(s % 100 == 0) dumper.Dump(); + if(s % 10 == 0) { + dumper.Dump(); + // if (psize > 1) dumper_ghost.Dump(); + } #endif //AKANTU_USE_IOHELPER - if(s % 10 == 0) std::cout << "passing step " << s << "/" << max_steps << std::endl; + if(prank == 0) + if(s % 10 == 0) std::cerr << "passing step " << s << "/" << max_steps << std::endl; } -#ifdef CHECK_STRESS - outfile.close(); -#endif // CHECK_STRESS + delete [] part; + delete model; akantu::finalize(); - return EXIT_SUCCESS; } diff --git a/test/test_solid_mechanics_model/test_solid_mechanics_model_bar_traction2d_parallel.cc b/test/test_solid_mechanics_model/test_solid_mechanics_model_line_parallel.cc similarity index 54% copy from test/test_solid_mechanics_model/test_solid_mechanics_model_bar_traction2d_parallel.cc copy to test/test_solid_mechanics_model/test_solid_mechanics_model_line_parallel.cc index f6d133faf..5102a0e02 100644 --- a/test/test_solid_mechanics_model/test_solid_mechanics_model_bar_traction2d_parallel.cc +++ b/test/test_solid_mechanics_model/test_solid_mechanics_model_line_parallel.cc @@ -1,228 +1,210 @@ /** * @file test_solid_mechanics_model.cc * @author Nicolas Richart * @date Tue Jul 27 14:34:13 2010 * * @brief test of the class SolidMechanicsModel * * @section LICENSE * * * */ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "static_communicator.hh" #include "communicator.hh" #include "mesh_partition_scotch.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER #define CHECK_STRESS - int main(int argc, char *argv[]) { - akantu::UInt spatial_dimension = 2; - akantu::UInt max_steps = 1000; + akantu::UInt spatial_dimension = 1; + akantu::ElementType type = akantu::_line_1; + akantu::UInt max_steps = 10000; akantu::Real time_factor = 0.2; akantu::initialize(&argc, &argv); - akantu::debug::setDebugLevel(akantu::dblDump); + // akantu::debug::setDebugLevel(akantu::dblDump); // akantu::Real epot, ekin; akantu::Mesh mesh(spatial_dimension); akantu::StaticCommunicator * comm = akantu::StaticCommunicator::getStaticCommunicator(); akantu::Int psize = comm->getNbProc(); akantu::Int prank = comm->whoAmI(); akantu::Communicator * communicator; if(prank == 0) { akantu::MeshIOMSH mesh_io; - mesh_io.read("bar.msh", mesh); - akantu::MeshPartition * partition = + mesh_io.read("line.msh", mesh); + akantu::MeshPartition * partition = new akantu::MeshPartitionScotch(mesh, spatial_dimension); partition->partitionate(psize); communicator = akantu::Communicator::createCommunicatorDistributeMesh(mesh, partition); delete partition; } else { communicator = akantu::Communicator::createCommunicatorDistributeMesh(mesh, NULL); } + // std::cout << mesh << std::endl; + akantu::SolidMechanicsModel * model = new akantu::SolidMechanicsModel(mesh); akantu::UInt nb_nodes = model->getFEM().getMesh().getNbNodes(); - akantu::UInt nb_element = model->getFEM().getMesh().getNbElement(akantu::_triangle_1); + akantu::UInt nb_element = model->getFEM().getMesh().getNbElement(type); /// model initialization model->initVectors(); /// set vectors to 0 memset(model->getForce().values, 0, spatial_dimension*nb_nodes*sizeof(akantu::Real)); memset(model->getVelocity().values, 0, spatial_dimension*nb_nodes*sizeof(akantu::Real)); memset(model->getAcceleration().values, 0, spatial_dimension*nb_nodes*sizeof(akantu::Real)); memset(model->getDisplacement().values, 0, spatial_dimension*nb_nodes*sizeof(akantu::Real)); model->readMaterials("material.dat"); model->initMaterials(); model->registerSynchronizer(*communicator); model->initModel(); std::cout << model->getMaterial(0) << std::endl; model->assembleMass(); #ifdef AKANTU_USE_IOHELPER /// set to 0 only for the first paraview dump memset(model->getResidual().values, 0, spatial_dimension*nb_nodes*sizeof(akantu::Real)); - memset(model->getMaterial(0).getStrain(akantu::_triangle_1).values, 0, + memset(model->getMaterial(0).getStrain(type).values, 0, spatial_dimension*spatial_dimension*nb_element*sizeof(akantu::Real)); - memset(model->getMaterial(0).getStress(akantu::_triangle_1).values, 0, + memset(model->getMaterial(0).getStress(type).values, 0, spatial_dimension*spatial_dimension*nb_element*sizeof(akantu::Real)); #endif //AKANTU_USE_IOHELPER /// boundary conditions - akantu::Real eps = 1e-16; for (akantu::UInt i = 0; i < nb_nodes; ++i) { - if(model->getFEM().getMesh().getNodes().values[spatial_dimension*i] >= 9) - model->getDisplacement().values[spatial_dimension*i] = (model->getFEM().getMesh().getNodes().values[spatial_dimension*i] - 9) / 100. ; - - if(model->getFEM().getMesh().getNodes().values[spatial_dimension*i] <= eps) - model->getBoundary().values[spatial_dimension*i] = true; + model->getDisplacement().values[spatial_dimension*i] = model->getFEM().getMesh().getNodes().values[i] / 100. ; - if(model->getFEM().getMesh().getNodes().values[spatial_dimension*i + 1] <= eps || - model->getFEM().getMesh().getNodes().values[spatial_dimension*i + 1] >= 1 - eps ) { - model->getBoundary().values[spatial_dimension*i + 1] = true; - } + if(model->getFEM().getMesh().getNodes().values[spatial_dimension*i] <= 1e-15) + model->getBoundary().values[i] = true; } + memset(model->getForce().values, 0, + spatial_dimension*nb_nodes*sizeof(akantu::Real)); + + // model->synchronizeBoundaries(); + + akantu::Real time_step = model->getStableTimeStep() * time_factor; std::cout << "Time Step = " << time_step << "s" << std::endl; model->setTimeStep(time_step); - // model->setTimeStep(3.54379e-07); #ifdef AKANTU_USE_IOHELPER DumperParaview dumper; dumper.SetMode(TEXT); + dumper.SetParallelContext(prank, psize); dumper.SetPoints(model->getFEM().getMesh().getNodes().values, - spatial_dimension, nb_nodes, "coordinates_par"); - dumper.SetConnectivity((int *)model->getFEM().getMesh().getConnectivity(akantu::_triangle_1).values, - TRIANGLE1, nb_element, C_MODE); + spatial_dimension, nb_nodes, "line_para"); + dumper.SetConnectivity((int *)model->getFEM().getMesh().getConnectivity(type).values, + LINE1, nb_element, C_MODE); dumper.AddNodeDataField(model->getDisplacement().values, spatial_dimension, "displacements"); dumper.AddNodeDataField(model->getVelocity().values, spatial_dimension, "velocity"); dumper.AddNodeDataField(model->getResidual().values, spatial_dimension, "force"); - dumper.AddElemDataField(model->getMaterial(0).getStrain(akantu::_triangle_1).values, + dumper.AddNodeDataField(model->getAcceleration().values, + spatial_dimension, "acceleration"); + dumper.AddElemDataField(model->getMaterial(0).getStrain(type).values, spatial_dimension*spatial_dimension, "strain"); - dumper.AddElemDataField(model->getMaterial(0).getStress(akantu::_triangle_1).values, + dumper.AddElemDataField(model->getMaterial(0).getStress(type).values, spatial_dimension*spatial_dimension, "stress"); - dumper.SetEmbeddedValue("displacements", 1); + double * part = new double[nb_element]; + for (unsigned int i = 0; i < nb_element; ++i) + part[i] = prank; + dumper.AddElemDataField(part, 1, "partitions"); + // dumper.SetEmbeddedValue("displacements", 1); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); -#endif //AKANTU_USE_IOHELPER + unsigned int nb_ghost_element = mesh.getNbGhostElement(type); + DumperParaview dumper_ghost; + dumper_ghost.SetMode(TEXT); + dumper_ghost.SetParallelContext(prank, psize); + dumper_ghost.SetPoints(mesh.getNodes().values, + spatial_dimension, nb_nodes, "line_para_ghost"); + dumper_ghost.SetConnectivity((int*) mesh.getGhostConnectivity(type).values, + LINE1, nb_ghost_element, C_MODE); + dumper_ghost.AddNodeDataField(model->getDisplacement().values, + spatial_dimension, "displacements"); + dumper_ghost.AddNodeDataField(model->getVelocity().values, + spatial_dimension, "velocity"); + dumper_ghost.AddNodeDataField(model->getResidual().values, + spatial_dimension, "force"); + dumper_ghost.AddNodeDataField(model->getAcceleration().values, + spatial_dimension, "acceleration"); + + dumper_ghost.SetPrefix("paraview/"); + dumper_ghost.Init(); + dumper_ghost.Dump(); +#endif //AKANTU_USE_IOHELPER -#ifdef CHECK_STRESS - std::ofstream outfile; - outfile.open("stress"); -#endif // CHECK_STRESS model->setPotentialEnergyFlagOn(); for(akantu::UInt s = 1; s <= max_steps; ++s) { model->explicitPred(); model->updateResidual(); model->updateAcceleration(); model->explicitCorr(); - // epot = model->getPotentialEnergy(); - // ekin = model->getKineticEnergy(); - // std::cout << s << " " << epot << " " << ekin << " " << epot + ekin - // << std::endl; - -#ifdef CHECK_STRESS - akantu::Real max_stress = std::numeric_limits::min(); - akantu::UInt max_el = 0; - akantu::Real * stress = model->getMaterial(0).getStress(akantu::_triangle_1).values; - for (akantu::UInt i = 0; i < nb_element; ++i) { - if(max_stress < stress[i*spatial_dimension*spatial_dimension]) { - max_stress = stress[i*spatial_dimension*spatial_dimension]; - max_el = i; - } - } + akantu::Real epot = model->getPotentialEnergy(); + akantu::Real ekin = model->getKineticEnergy(); - akantu::Real * coord = model->getFEM().getMesh().getNodes().values; - akantu::Real * disp_val = model->getDisplacement().values; - akantu::UInt * conn = model->getFEM().getMesh().getConnectivity(akantu::_triangle_1).values; - akantu::UInt nb_nodes_per_element = model->getFEM().getMesh().getNbNodesPerElement(akantu::_triangle_1); - akantu::Real * coords = new akantu::Real[spatial_dimension]; - akantu::Real min_x = std::numeric_limits::max(); - akantu::Real max_x = std::numeric_limits::min(); - - akantu::Real stress_range = 5e7; - for (akantu::UInt el = 0; el < nb_element; ++el) { - if(stress[el*spatial_dimension*spatial_dimension] > max_stress - stress_range) { - akantu::UInt el_offset = el * nb_nodes_per_element; - memset(coords, 0, spatial_dimension*sizeof(akantu::Real)); - for (akantu::UInt n = 0; n < nb_nodes_per_element; ++n) { - for (akantu::UInt i = 0; i < spatial_dimension; ++i) { - akantu::UInt node = conn[el_offset + n] * spatial_dimension; - coords[i] += (coord[node + i] + disp_val[node + i]) - / ((akantu::Real) nb_nodes_per_element); - } - } - min_x = min_x < coords[0] ? min_x : coords[0]; - max_x = max_x > coords[0] ? max_x : coords[0]; - } + if(prank == 0) { + std::cout << s << " " << epot << " " << ekin << " " << epot + ekin + << std::endl; } - - - outfile << s << " " << .5 * (min_x + max_x) << " " << min_x << " " << max_x << " " << max_x - min_x << " " << max_stress << std::endl; - - delete [] coords; -#endif // CHECK_STRESS - - #ifdef AKANTU_USE_IOHELPER - if(s % 100 == 0) dumper.Dump(); + // if(s % 100 == 0) + dumper.Dump(); + dumper_ghost.Dump(); #endif //AKANTU_USE_IOHELPER - if(s % 10 == 0) std::cout << "passing step " << s << "/" << max_steps << std::endl; + if(s % 10 == 0) std::cerr << "passing step " << s << "/" << max_steps << std::endl; } -#ifdef CHECK_STRESS - outfile.close(); -#endif // CHECK_STRESS + delete [] part; akantu::finalize(); return EXIT_SUCCESS; }