diff --git a/VERSION b/VERSION index 60a8ddb66..858ceb575 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -1.0.0.5964 +1.0.0.5973 diff --git a/cmake/AkantuBuildTreeSettings.cmake.in b/cmake/AkantuBuildTreeSettings.cmake.in index c91764a06..ab43df496 100644 --- a/cmake/AkantuBuildTreeSettings.cmake.in +++ b/cmake/AkantuBuildTreeSettings.cmake.in @@ -1,4 +1,9 @@ set(AKANTU_INCLUDE_DIR "@AKANTU_INCLUDE_DIRS@" ) + +list(APPEND AKANTU_OPTION_LIST IOHELPER) +set(AKANTU_IOHELPER_INCLUDE_DIR @CMAKE_SOURCE_DIR@/third-party/iohelper/src) +#set(AKANTU_IOHELPER_LIBRARIES @CMAKE_BINARY_DIR@/third-party/iohelper/libiohelper.so) + set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "@CMAKE_SOURCE_DIR@/cmake") \ No newline at end of file diff --git a/cmake/AkantuInstall.cmake b/cmake/AkantuInstall.cmake index 648f57c3b..977f136c9 100644 --- a/cmake/AkantuInstall.cmake +++ b/cmake/AkantuInstall.cmake @@ -1,119 +1,119 @@ #=============================================================================== # @file AkantuInstall.cmake # @author Nicolas Richart # @date Wed Oct 17 2012 # # @brief Create the files that allows users to link with Akantu in an other # cmake project # # @section LICENSE # # Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # #=============================================================================== #=============================================================================== # Config gen for external packages #=============================================================================== configure_file(cmake/AkantuBuildTreeSettings.cmake.in "${CMAKE_BINARY_DIR}/AkantuBuildTreeSettings.cmake" @ONLY) file(WRITE "${CMAKE_BINARY_DIR}/AkantuConfigInclude.cmake" " #=============================================================================== # @file AkantuConfigInclude.cmake # @author Nicolas Richart # @date Fri Jun 11 09:46:59 2010 # # @section LICENSE # # Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== ") foreach(_option ${PACKAGE_SYSTEM_PACKAGES_NAMES_LIST_ALL}) - list(FIND AKANTU_OPTION_LIST ${_option} _index) + list(FIND AKANTU_OPTION_LIST ${_option_name} _index) if (_index EQUAL -1) if(NOT "${_option}" STREQUAL "CORE") if(NOT AKANTU_${_option}) set(AKANTU_${_option} OFF) endif() file(APPEND "${CMAKE_BINARY_DIR}/AkantuConfigInclude.cmake" " set(AKANTU_HAS_${_option} ${AKANTU_${_option}})") endif() endif() endforeach() file(APPEND "${CMAKE_BINARY_DIR}/AkantuConfigInclude.cmake" " set(AKANTU_HAS_PARTITIONER ${AKANTU_PARTITIONER}) set(AKANTU_HAS_SOLVER ${AKANTU_SOLVER}) - -set(AKANTU_OPTION_LIST ${AKANTU_OPTION_LIST}) ") foreach(_option ${AKANTU_OPTION_LIST}) + package_pkg_name(${_option} _pkg_name) file(APPEND "${CMAKE_BINARY_DIR}/AkantuConfigInclude.cmake" " +list(APPEND AKANTU_OPTION_LIST ${_option}) set(AKANTU_USE_${_option} ${AKANTU_${_option}})") - if(${AKANTU_${_option}_LIBRARIES}) + if(${_pkg_name}_LIBRARIES) file(APPEND "${CMAKE_BINARY_DIR}/AkantuConfigInclude.cmake" " -set(AKANTU_${_option}_LIBRARIES ${AKANTU_${_option}_LIBRARIES}) -set(AKANTU_${_option}_INCLUDE_DIR ${AKANTU_${_option}_INCLUDE_DIR}) +set(${_pkg_name}_LIBRARIES ${${_pkg_name}_LIBRARIES}) +set(${_pkg_name}_INCLUDE_DIR ${${_pkg_name}_INCLUDE_DIR}) ") endif() endforeach() # Create the AkantuConfig.cmake and AkantuConfigVersion files get_filename_component(CONF_REL_INCLUDE_DIR "${CMAKE_INSTALL_PREFIX}" ABSOLUTE) configure_file(cmake/AkantuConfig.cmake.in "${CMAKE_BINARY_DIR}/AkantuConfig.cmake" @ONLY) configure_file(cmake/AkantuConfigVersion.cmake.in "${CMAKE_BINARY_DIR}/AkantuConfigVersion.cmake" @ONLY) configure_file(cmake/AkantuUse.cmake "${CMAKE_BINARY_DIR}/AkantuUse.cmake" COPYONLY) # Install the export set for use with the install-tree install(FILES ${CMAKE_BINARY_DIR}/AkantuConfig.cmake ${CMAKE_BINARY_DIR}/AkantuConfigInclude.cmake ${CMAKE_BINARY_DIR}/AkantuConfigVersion.cmake ${CMAKE_SOURCE_DIR}/cmake/AkantuUse.cmake DESTINATION lib/akantu COMPONENT dev) install(FILES ${CMAKE_SOURCE_DIR}/cmake/FindIOHelper.cmake ${CMAKE_SOURCE_DIR}/cmake/FindQVIEW.cmake ${CMAKE_SOURCE_DIR}/cmake/FindMumps.cmake ${CMAKE_SOURCE_DIR}/cmake/FindScotch.cmake ${CMAKE_SOURCE_DIR}/cmake/FindGMSH.cmake DESTINATION lib/akantu/cmake COMPONENT dev) diff --git a/cmake/AkantuMacros.cmake b/cmake/AkantuMacros.cmake index 6ac8b5192..facb60662 100644 --- a/cmake/AkantuMacros.cmake +++ b/cmake/AkantuMacros.cmake @@ -1,154 +1,154 @@ #=============================================================================== # @file AkantuMacros.cmake # @author Nicolas Richart # @author Guillaume Anciaux # @date Wed Feb 9 10:59:42 2011 # # @brief Set of macros used by akantu cmake files # # @section LICENSE # # Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # #=============================================================================== -set(AKANTU_CMAKE_DEBUG FALSE) +set(AKANTU_CMAKE_DEBUG TRUE) macro(akantu_message) if(AKANTU_CMAKE_DEBUG) message(${ARGN}) endif() endmacro() #=============================================================================== macro(add_cxx_flags flag) if(NOT CMAKE_CXX_FLAGS MATCHES "${flag}") set(CMAKE_CXX_FLAGS "${flag} ${CMAKE_CXX_FLAGS}" CACHE STRING "Flags used by the compiler during all build types." FORCE) endif() endmacro() #=============================================================================== macro(remove_cxx_flags flag) string(REPLACE "${flag} " "" CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS}") set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS}" CACHE STRING "Flags used by the compiler during all build types." FORCE) endmacro() #=============================================================================== macro(check_for_isnan result) include(CheckFunctionExists) check_function_exists(std::isnan HAVE_STD_ISNAN) if(HAVE_STD_ISNAN) set(result "std::isnan(x)") else() check_function_exists(isnan HAVE_ISNAN) if(HAVE_ISNAN) set(result "(::isnan(x))") else() check_function_exists(_isnan HAVE_ISNAN_MATH_H) if(HAVE_ISNAN_MATH_H) set(result "(_isnan(x))") else() set(result (x == std::numeric_limits::quiet_NAN())) endif() endif() endif() endmacro() #=============================================================================== macro(copy_files target_depend) foreach(_file ${ARGN}) set(_target ${CMAKE_CURRENT_BINARY_DIR}/${_file}) set(_source ${CMAKE_CURRENT_SOURCE_DIR}/${_file}) add_custom_command( OUTPUT ${_target} COMMAND ${CMAKE_COMMAND} -E copy_if_different ${_source} ${_target} DEPENDS ${_source} ) set(_target_name ${target_depend}_${_file}) add_custom_target(${_target_name} DEPENDS ${_target}) add_dependencies(${target_depend} ${_target_name}) endforeach() endmacro() #=============================================================================== macro(add_akantu_definitions) foreach(_definition ${AKANTU_DEFINITIONS}) add_definitions(-D${_definition}) endforeach() endmacro() #=============================================================================== if(__CMAKE_PARSE_ARGUMENTS_INCLUDED) return() endif() set(__CMAKE_PARSE_ARGUMENTS_INCLUDED TRUE) function(CMAKE_PARSE_ARGUMENTS prefix _optionNames _singleArgNames _multiArgNames) # first set all result variables to empty/FALSE foreach(arg_name ${_singleArgNames} ${_multiArgNames}) set(${prefix}_${arg_name}) endforeach(arg_name) foreach(option ${_optionNames}) set(${prefix}_${option} FALSE) endforeach(option) set(${prefix}_UNPARSED_ARGUMENTS) set(insideValues FALSE) set(currentArgName) # now iterate over all arguments and fill the result variables foreach(currentArg ${ARGN}) list(FIND _optionNames "${currentArg}" optionIndex) # ... then this marks the end of the arguments belonging to this keyword list(FIND _singleArgNames "${currentArg}" singleArgIndex) # ... then this marks the end of the arguments belonging to this keyword list(FIND _multiArgNames "${currentArg}" multiArgIndex) # ... then this marks the end of the arguments belonging to this keyword if(${optionIndex} EQUAL -1 AND ${singleArgIndex} EQUAL -1 AND ${multiArgIndex} EQUAL -1) if(insideValues) if("${insideValues}" STREQUAL "SINGLE") set(${prefix}_${currentArgName} ${currentArg}) set(insideValues FALSE) elseif("${insideValues}" STREQUAL "MULTI") list(APPEND ${prefix}_${currentArgName} ${currentArg}) endif() else(insideValues) list(APPEND ${prefix}_UNPARSED_ARGUMENTS ${currentArg}) endif(insideValues) else() if(NOT ${optionIndex} EQUAL -1) set(${prefix}_${currentArg} TRUE) set(insideValues FALSE) elseif(NOT ${singleArgIndex} EQUAL -1) set(currentArgName ${currentArg}) set(${prefix}_${currentArgName}) set(insideValues "SINGLE") elseif(NOT ${multiArgIndex} EQUAL -1) set(currentArgName ${currentArg}) set(${prefix}_${currentArgName}) set(insideValues "MULTI") endif() endif() endforeach(currentArg) # propagate the result variables to the caller: foreach(arg_name ${_singleArgNames} ${_multiArgNames} ${_optionNames}) set(${prefix}_${arg_name} ${${prefix}_${arg_name}} PARENT_SCOPE) endforeach(arg_name) set(${prefix}_UNPARSED_ARGUMENTS ${${prefix}_UNPARSED_ARGUMENTS} PARENT_SCOPE) endfunction(CMAKE_PARSE_ARGUMENTS _options _singleArgs _multiArgs) diff --git a/cmake/AkantuPackageSystem.cmake b/cmake/AkantuPackageSystem.cmake index fe0adee07..f87163c84 100644 --- a/cmake/AkantuPackageSystem.cmake +++ b/cmake/AkantuPackageSystem.cmake @@ -1,279 +1,299 @@ #=============================================================================== # @file AkantuMacros.cmake # @author Nicolas Richart # @author Guillaume Anciaux # @date Wed Feb 9 10:59:42 2011 # # @brief Set of macros used by akantu to handle the package system # # @section LICENSE # # Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # #=============================================================================== #=============================================================================== # Package Management #=============================================================================== set(PACKAGE_SYSTEM_PKG_PREFIX AKANTU_) set(PACKAGE_SYSTEM_OPT_PREFIX AKANTU_USE_) set(PACKAGE_SYSTEM_DESC_PREFIX AKANTU_DESC_) macro(package_pkg_name PKG PKG_NAME) string(TOUPPER ${PKG} _u_package) set(${PKG_NAME} ${PACKAGE_SYSTEM_PKG_PREFIX}${_u_package}) endmacro() macro(package_opt_name PKG OPT_NAME) string(TOUPPER ${PKG} _u_package) set(${OPT_NAME} ${PACKAGE_SYSTEM_OPT_PREFIX}${_u_package}) endmacro() macro(package_desc_name PKG DESC_NAME) string(TOUPPER ${PKG} _u_package) set(${DESC_NAME} ${PACKAGE_SYSTEM_DESC_PREFIX}${_u_package}) endmacro() #=============================================================================== macro(add_all_packages package_dir) akantu_message("add_all_packages: PKG DIR : ${package_dir}") file(GLOB _akantu_package_list "${package_dir}/*.cmake") foreach(_pkg ${_akantu_package_list}) get_filename_component(_basename ${_pkg} NAME) string(REGEX REPLACE "\\.cmake" "" _option_name ${_basename}) - string(TOUPPER "${_pkg_suffix}" _pkg_suffix) + string(TOUPPER "${_option_name}" _option_name) list(APPEND PACKAGE_SYSTEM_PACKAGES_NAMES_LIST_ALL ${_option_name}) endforeach() akantu_message("add_all_packages: PKG LIST : "${PACKAGE_SYSTEM_PACKAGES_NAMES_LIST_ALL}) foreach(_pkg ${PACKAGE_SYSTEM_PACKAGES_NAMES_LIST_ALL}) akantu_message("add_all_packages: including ${_pkg}") + string(TOLOWER "${_pkg}" _pkg) include(${package_dir}/${_pkg}.cmake) package_pkg_name(${_pkg} _package_name) if (${_package_name}) list(APPEND PACKAGE_SYSTEM_PACKAGES_ON ${_pkg}) else (${_package_name}) list(APPEND PACKAGE_SYSTEM_PACKAGES_OFF ${_pkg}) endif() foreach(_file ${${_pkg_name}_FILES}) list(APPEND _release_all_file ${_file}) endforeach() endforeach() akantu_message("add_all_packages: ON PKG : ${PACKAGE_SYSTEM_PACKAGES_ON}") akantu_message("add_all_packages: ALL FILE LIST : ${_release_all_file}") #check if there are some file in the release that are not registered in a package file(GLOB_RECURSE _all_files RELATIVE ${CMAKE_SOURCE_DIR}/src "*.cc" "*.hh") foreach(_file ${all_files}) if(NOT ${_file} MATCHES "test.*") list(FIND _release_all_file ${_file} _index) if (_index EQUAL -1) message("The file ${file} is not registered in any package.") message("Please append the file in one of the files within directory ${CMAKE_SOURCE_DIR}/packages") endif() endif() endforeach() #check if there are some file in the package list that are not on the current directory foreach(file ${all_files}) list(FIND RELEASE_ALL_FILE ${file} index) if (index EQUAL -1) message("file ${file} is not registered in any package.") message("Please append the file in one of the files within directory ${AKANTU_SOURCE_DIR}/packages") endif() endforeach() foreach(_file ${_release_all_file}) list(FIND _all_files ${_file} _index) if (_index EQUAL -1) message("The file ${_file} is registered in packages but is not present in the source directory.") endif() endforeach() #construct list of files for unactivated packages foreach(_pkg ${PACKAGE_SYSTEM_PACKAGES_OFF}) package_pkg_name(${_pkg} _pkg_name) foreach(_file ${${_pkg_name}_FILES}) # string(REGEX REPLACE "\\/" "\\\\\\\\/" __file ${_file}) akantu_message("add_all_packages: ${_file} ${__file}") # list(APPEND _exclude_source_file "/${__file}/") list(APPEND AKANTU_EXCLUDE_SOURCE_FILE ${_file}) endforeach() endforeach() #check dependencies foreach(_pkg ${PACKAGE_SYSTEM_PACKAGES_OFF}) # differentiate the file types akantu_message("add_all_packages: DEPENDS PKG : ${_pkg}") akantu_message("add_all_packages: DEPENDS LST : ${${_pkg}_DEPENDS}") package_pkg_name(${_pkg} _pkg_name) if (NOT "${${_pkg_name}_DEB_DEPEND}" STREQUAL "") set(deb_depend "${deb_depend}, ${${_pkg}_DEB_DEPEND}") endif() endforeach() set(PACKAGE_SYSTEM_DEBIAN_PACKAGE_DEPENDS "${deb_depend}") endmacro() #=============================================================================== macro(generate_source_list_from_packages source_dir source_files headers_files include_dirs) set(deb_depend "libc6") akantu_message("generate_source_list_from_packages: SRC DIR : ${source_dir}") foreach(_pkg ${PACKAGE_SYSTEM_PACKAGES_ON}) # differentiate the file types package_pkg_name(${_pkg} _package_name) foreach(_file ${${_package_name}_FILES}) if(${_file} MATCHES ".*inline.*\\.cc") # list(APPEND ${_package_name}_inlines ${_file}) list(APPEND ${_package_name}_headers ${_file}) elseif(${_file} MATCHES ".*\\.hh") list(APPEND ${_package_name}_headers ${_file}) else() list(APPEND ${_package_name}_srcs ${_file}) endif() endforeach() # generates the include directory variable foreach(_file ${${_package_name}_headers}) get_filename_component(_absolute_name ${_file} ABSOLUTE) get_filename_component(_include_dir ${_absolute_name} PATH) list(APPEND ${_package_name}_include_dirs ${_include_dir}) list(REMOVE_DUPLICATES ${_package_name}_include_dirs) endforeach() # generate global lists for akantu to know what to build list(APPEND ${source_files} ${${_package_name}_srcs}) list(APPEND ${headers_files} ${${_package_name}_headers}) list(APPEND ${include_dirs} ${${_package_name}_include_dirs}) akantu_message("generate_source_list_from_packages: PKG ${_package_name} SRCS : ${${_package_name}_srcs}") akantu_message("generate_source_list_from_packages: PKG ${_package_name} HRDS : ${${_package_name}_headers}") akantu_message("generate_source_list_from_packages: PKG ${_package_name} INCS : ${${_package_name}_include_dirs}") endforeach() akantu_message("generate_source_list_from_packages: SRCS : ${${source_files}}") akantu_message("generate_source_list_from_packages: HRDS : ${${headers_files}}") akantu_message("generate_source_list_from_packages: INCS : ${${include_dirs}}") endmacro() #=============================================================================== # macro to include optional packages macro(add_optional_package PACKAGE DESC DEFAULT) cmake_parse_arguments(_opt_pkg "" "LANGUAGE" "DEPENDS;PREFIX" ${ARGN}) package_pkg_name (${PACKAGE} _pkg_name) package_opt_name (${PACKAGE} _option_name) package_desc_name(${PACKAGE} _desc_name) akantu_message("add_optional_package: Registering ${PACKAGE} ${DESC} -> ${_option_name}") if(_opt_pkg_PREFIX) set(_package_prefix ${_opt_pkg_PREFIX}) else() string(TOUPPER ${PACKAGE} _u_package) set(_package_prefix ${_u_package}) endif() set(${_desc_name} ${DESC}) option(${_option_name} ${DESC} ${DEFAULT}) if(${_option_name}) if(_opt_pkg_LANGUAGE) foreach(_language ${_opt_pkg_LANGUAGE}) akantu_message("add_optional_package: Package ${PACKAGE} asked for language ${_language}") enable_language(${_language}) endforeach() endif() foreach(_dep ${_opt_pkg_DEPENDS}) add_package_dependecies(${PACKAGE} ${_dep}) endforeach() find_package(${PACKAGE} REQUIRED) foreach(_prefix ${_package_prefix}) if(${_prefix}_FOUND) list(APPEND AKANTU_DEFINITIONS ${_option_name}) if(DEFINED ${_prefix}_INCLUDE_DIR) list(APPEND AKANTU_EXTERNAL_LIB_INCLUDE_DIR ${${_prefix}_INCLUDE_DIR}) set(${_pkg_name}_INCLUDE_DIR ${${_prefix}_INCLUDE_DIR}) else() list(APPEND AKANTU_EXTERNAL_LIB_INCLUDE_DIR ${${_prefix}_INCLUDE_PATH}) set(${_pkg_name}_INCLUDE_DIR ${${_prefix}_INCLUDE_PATH}) endif() list(APPEND AKANTU_EXTERNAL_LIBRARIES ${${_prefix}_LIBRARIES}) set(${_pkg_name}_LIBRARIES ${${_prefix}_LIBRARIES}) set(${_pkg_name} ON) + string(TOUPPER ${PACKAGE} _u_package) list(APPEND AKANTU_OPTION_LIST ${_u_package}) - akantu_message("add_optional_package: Package ${PACKAGE} found!") - akantu_message("add_optional_package: Package ${PACKAGE} includes : ${_pkg_name}_INCLUDE_DIR}") - akantu_message("add_optional_package: Package ${PACKAGE} libraries: ${_pkg_name}_LIBRARIES}") + akantu_message("add_optional_package: Package ${PACKAGE} includes : ${${_pkg_name}_INCLUDE_DIR}") + akantu_message("add_optional_package: Package ${PACKAGE} libraries: ${${_pkg_name}_LIBRARIES}") + akantu_message("add_optional_package: option list: ${AKANTU_OPTION_LIST}") else(${_prefix}_FOUND) akantu_message("add_optional_package: Package ${PACKAGE} not found!") set(${_pkg_name} OFF) endif(${_prefix}_FOUND) endforeach() endif(${_option_name}) endmacro() #=============================================================================== macro(add_package_dependecies PKG DEP) - package_opt_name (${PKG} _opt_name) + package_pkg_name (${PKG} _opt_name) package_opt_name (${DEP} _dep_name) package_desc_name(${DEP} _var_dep_desc) akantu_message("add_package_dependecies: add dependence between ${_opt_name} and ${_dep_name}") - set(_dep_desc ${_var_dep_desc}) + akantu_message("add_package_dependecies: ON dependencies of ${_dep_name} are: ${${_dep_name}_DEPS}") + akantu_message("add_package_dependecies: saved value for ${_dep_name} is: ${${_dep_name}_OLD}") if(${_opt_name}) - akantu_message("add_package_dependecies: Save dep state ${_dep_name}:${${_dep_name}}") - set(${_dep_name}_OLD ${_dep_name} CACHE INTERNAL "${_dep_desc}" FORCE) + if("${${_dep_name}_DEPS}" STREQUAL "") + akantu_message("add_package_dependecies: Save dep state ${_dep_name}:${${_dep_name}}") + set(${_dep_name}_OLD ${${_dep_name}} CACHE INTERNAL "${_dep_desc}" FORCE) + endif() + + akantu_message("add_package_dependecies: force value to ON ${_dep_name}") set(${_dep_name} ON CACHE BOOL "${_dep_desc}" FORCE) + + list(FIND ${_dep_name}_DEPS ${_opt_name} pos) + if(pos EQUAL -1) + list(APPEND ${_dep_name}_DEPS ${_opt_name}) + set(${_dep_name}_DEPS ${${_dep_name}_DEPS} CACHE INTERNAL "Dependencies ON with package ${_dep_name}" FORCE) + endif() else() - if(DEFINED ${_dep_name}_OLD) - akantu_message("add_package_dependecies: Restore state ${_dep_name}:${${_dep_name}}") - set(${_dep_name} ${_dep_name}_OLD CACHE BOOL "${_dep_desc}" FORCE) + list(LENGTH ${_dep_name}_DEPS len) + list(FIND ${_dep_name}_DEPS ${_opt_name} pos) + if((len EQUAL 1) AND (NOT pos EQUAL -1)) + akantu_message("add_package_dependecies: Restore state ${_dep_name}:${${_dep_name}} (${pos})") + set(${_dep_name} ${${_dep_name}_OLD} CACHE BOOL "${_dep_desc}" FORCE) unset(${_dep_name}_OLD CACHE) - endif(DEFINED ${_dep_name}_OLD) + endif() + + if(NOT pos EQUAL -1) + list(REMOVE_AT ${_dep_name}_DEPS ${pos}) + set(${_dep_name}_DEPS ${${_dep_name}_DEPS} CACHE INTERNAL "Nb dependencies with package ${_dep_name}" FORCE) + endif() endif() endmacro() #=============================================================================== # macro to add meta packages macro(add_meta_package PKG DESC DEFAULT) akantu_message("add_meta_package: register meta option ${PKG} ${DESC} ${DEFAULT}") package_pkg_name (${PKG} _pkg_name) package_desc_name(${PKG} _desc_name) set(${_desc_name} ${DESC}) option(${_pkg_name} ${DESC} ${DEFAULT}) foreach(_dep ${ARGN}) package_opt_name (${_dep} _dep_name) mark_as_advanced(${_dep_name}) add_package_dependecies(${PKG} ${_dep}) endforeach() endmacro() \ No newline at end of file diff --git a/packages/blas_lapack.cmake b/packages/blas.cmake similarity index 62% rename from packages/blas_lapack.cmake rename to packages/blas.cmake index aafd658f1..f6f01b3af 100644 --- a/packages/blas_lapack.cmake +++ b/packages/blas.cmake @@ -1,6 +1,5 @@ add_optional_package(BLAS "Use BLAS for arithmetic operations" OFF LANGUAGE Fortran) if(BLAS_mkl_core_LIBRARY) set(AKANTU_USE_BLAS_MKL) endif() -add_optional_package(LAPACK "Use LAPACK for arithmetic operations" OFF LANGUAGE Fortran) diff --git a/packages/iohelper.cmake b/packages/iohelper.cmake index 047bf45ba..9170b9941 100644 --- a/packages/iohelper.cmake +++ b/packages/iohelper.cmake @@ -1,13 +1,19 @@ #add_optional_package(IOHelper "Add IOHelper support in akantu" ON) option(AKANTU_USE_IOHELPER "Add IOHelper support in akantu" ON) mark_as_advanced(AKANTU_USE_IOHELPER) if(AKANTU_USE_IOHELPER) set(IOHELPER_TARGETS_EXPORT AkantuLibraryDepends) add_subdirectory(third-party/iohelper) + list(APPEND AKANTU_EXTERNAL_LIBRARIES iohelper) - mark_as_advanced(IOHELPER_TESTS) -endif() + list(APPEND AKANTU_EXTERNAL_LIB_INCLUDE_DIR ${AKANTU_BUILD_IOHELPER_INCLUDE_DIR}) + list(APPEND AKANTU_EXPORT_LIST iohelper) + mark_as_advanced(IOHELPER_TESTS) + set(AKANTU_IOHELPER ON) +else() + set(AKANTU_IOHELPER OFF) +endif() diff --git a/packages/lapack.cmake b/packages/lapack.cmake new file mode 100644 index 000000000..2d545a32b --- /dev/null +++ b/packages/lapack.cmake @@ -0,0 +1 @@ +add_optional_package(LAPACK "Use LAPACK for arithmetic operations" OFF LANGUAGE Fortran) \ No newline at end of file diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index a2d745958..02248fa4f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,97 +1,99 @@ #=============================================================================== # @file CMakeLists.txt # @author Nicolas Richart # @date Mon Nov 28 15:11:51 2011 # # @brief CMake file for the library # # @section LICENSE # # Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # #=============================================================================== #=============================================================================== # Package Management #=============================================================================== generate_source_list_from_packages(${CMAKE_CURRENT_BINARY_DIR} AKANTU_LIBRARY_SRCS AKANTU_LIBRARY_HDRS AKANTU_LIBRARY_INCLUDE_DIRS) #=========================================================================== # header precompilation #=========================================================================== set(AKANTU_COMMON_HDR_TO_PRECOMPILE common/aka_vector.hh common/aka_math.hh common/aka_types.hh fem/element_class.hh ) set(AKANTU_PRECOMPILE_HDR_INCLUDE_DIRS ${CMAKE_CURRENT_BINARY_DIR}/common/ ${CMAKE_CURRENT_BINARY_DIR}/fem/ ) list(APPEND AKANTU_INCLUDE_DIRS ${CMAKE_CURRENT_BINARY_DIR} ${AKANTU_LIBRARY_INCLUDE_DIRS} ${AKANTU_PRECOMPILE_HDR_INCLUDE_DIRS}) set(AKANTU_INCLUDE_DIRS ${AKANTU_INCLUDE_DIRS} PARENT_SCOPE) include_directories(${AKANTU_INCLUDE_DIRS} ${AKANTU_EXTERNAL_LIB_INCLUDE_DIR}) if(CMAKE_COMPILER_IS_GNUCXX) include(${CMAKE_SOURCE_DIR}/cmake/PCHgcc.cmake) foreach(_header ${AKANTU_COMMON_HDR_TO_PRECOMPILE}) add_pch_rule(${_header} AKANTU_LIBRARY_SRCS) endforeach() elseif(CMAKE_COMPILER_IS_GNUCXX) endif() #=============================================================================== # Library generation #=============================================================================== configure_file(common/aka_config.hh.in "${CMAKE_CURRENT_BINARY_DIR}/aka_config.hh" @ONLY) list(APPEND AKANTU_LIBRARY_HDRS ${CMAKE_CURRENT_BINARY_DIR}/aka_config.hh) add_library(akantu ${AKANTU_LIBRARY_SRCS}) target_link_libraries(akantu ${AKANTU_EXTERNAL_LIBRARIES}) set_target_properties(akantu PROPERTIES ${AKANTU_LIBRARY_PROPERTIES}) set_target_properties(akantu PROPERTIES PUBLIC_HEADER "${AKANTU_LIBRARY_HDRS}") +list(APPEND AKANTU_EXPORT_LIST akantu) + install(TARGETS akantu EXPORT AkantuLibraryDepends LIBRARY DESTINATION lib COMPONENT shlib PUBLIC_HEADER DESTINATION include/akantu/ COMPONENT dev ) install(EXPORT AkantuLibraryDepends DESTINATION lib/akantu COMPONENT dev) #Export for build tree -#export(TARGETS akantu -# FILE "${CMAKE_BINARY_DIR}/AkantuLibraryDepends.cmake") -#export(PACKAGE Akantu) +export(TARGETS ${AKANTU_EXPORT_LIST} + FILE "${CMAKE_BINARY_DIR}/AkantuLibraryDepends.cmake") +export(PACKAGE Akantu) diff --git a/src/common/aka_math_tmpl.hh b/src/common/aka_math_tmpl.hh index 24ee1c3ed..4f2b7067c 100644 --- a/src/common/aka_math_tmpl.hh +++ b/src/common/aka_math_tmpl.hh @@ -1,662 +1,662 @@ /** * @file aka_math_inline_impl.cc * @author Nicolas Richart * @date Wed Jul 28 13:20:35 2010 * * @brief Implementation of the inline functions of the math toolkit * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ __END_AKANTU__ #include #include #ifdef AKANTU_USE_BLAS # ifndef AKANTU_USE_BLAS_MKL # include # else // AKANTU_USE_BLAS_MKL # include # endif //AKANTU_USE_BLAS_MKL #endif //AKANTU_USE_BLAS #ifdef AKANTU_USE_LAPACK extern "C" { void dgeev_(char* jobvl, char* jobvr, int* n, double* a, int* lda, double* wr, double* wi, double* vl, int* ldvl, double* vr, int* ldvr, double* work, int* lwork, int* info); // LU decomposition of a general matrix void dgetrf_(int* m, int *n, double* a, int* lda, int* ipiv, int* info); // generate inverse of a matrix given its LU decomposition void dgetri_(int* n, double* a, int* lda, int* ipiv, double* work, int* lwork, int* info); } #endif __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ inline void Math::matrix_vector(UInt m, UInt n, const Real * A, const Real * x, Real * y, Real alpha) { #ifdef AKANTU_USE_BLAS /// y = alpha*op(A)*x + beta*y cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n, alpha, A, n, x, 1, 0, y, 1); #else memset(y, 0, m*sizeof(Real)); for (UInt i = 0; i < m; ++i) { UInt A_i = i * n; for (UInt j = 0; j < n; ++j) { y[i] += A[A_i + j] * x[j]; } y[i] *= alpha; } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrixt_vector(UInt m, UInt n, const Real * A, const Real * x, Real * y, Real alpha) { #ifdef AKANTU_USE_BLAS /// y = alpha*op(A)*x + beta*y cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n, alpha, A, m, x, 1, 0, y, 1); #else memset(y, 0, m*sizeof(Real)); for (UInt i = 0; i < m; ++i) { for (UInt j = 0; j < n; ++j) { y[i] += A[i + j * m] * x[j]; } y[i] *= alpha; } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrix_matrix(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C, Real alpha) { #ifdef AKANTU_USE_BLAS /// C := alpha*op(A)*op(B) + beta*C cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, A, k, B, n, 0, C, n); #else memset(C, 0, m*n*sizeof(Real)); for (UInt j = 0; j < n; ++j) { for (UInt i = 0; i < m; ++i) { UInt A_i = i * k; UInt C_i = i * n; for (UInt l = 0; l < k; ++l) { UInt B_l = l * n; C[C_i + j] += A[A_i + l] * B[B_l + j]; } C[C_i + j] *= alpha; } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrixt_matrix(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C, Real alpha) { #ifdef AKANTU_USE_BLAS /// C := alpha*op(A)*op(B) + beta*C cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans, m, n, k, alpha, A, m, B, n, 0, C, n); #else memset(C, 0, m*n*sizeof(Real)); for (UInt i = 0; i < m; ++i) { // UInt A_i = i * k; UInt C_i = i * n; for (UInt j = 0; j < n; ++j) { for (UInt l = 0; l < k; ++l) { C[C_i + j] += A[l * m + i] * B[l * n + j]; } C[C_i + j] *= alpha; } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrix_matrixt(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C, Real alpha) { #ifdef AKANTU_USE_BLAS /// C := alpha*op(A)*op(B) + beta*C cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, m, n, k, alpha, A, k, B, k, 0, C, n); #else memset(C, 0, m*n*sizeof(Real)); for (UInt i = 0; i < m; ++i) { UInt A_i = i * k; UInt C_i = i * n; for (UInt j = 0; j < n; ++j) { UInt B_j = j * k; for (UInt l = 0; l < k; ++l) { C[C_i + j] += A[A_i + l] * B[B_j + l]; } C[C_i + j] *= alpha; } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrixt_matrixt(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C, Real alpha) { #ifdef AKANTU_USE_BLAS /// C := alpha*op(A)*op(B) + beta*C cblas_dgemm(CblasRowMajor, CblasTrans, CblasTrans, m, n, k, alpha, A, m, B, k, 0, C, n); #else memset(C, 0, m * n * sizeof(Real)); for (UInt i = 0; i < m; ++i) { UInt C_i = i * n; for (UInt j = 0; j < n; ++j) { UInt B_j = j * k; for (UInt l = 0; l < k; ++l) { C[C_i + j] += A[l * m + i] * B[B_j + l]; } C[C_i + j] *= alpha; } } #endif } /* -------------------------------------------------------------------------- */ inline Real Math::vectorDot(const Real * v1, const Real * v2, UInt n) { #ifdef AKANTU_USE_BLAS /// d := v1 . v2 - cblas_ddot(n, v1, 1, v2, 1); + Real d = cblas_ddot(n, v1, 1, v2, 1); #else Real d = 0; for (UInt i = 0; i < n; ++i) { d += v1[i] * v2[i]; } #endif return d; } /* -------------------------------------------------------------------------- */ template inline void Math::matMul(UInt m, UInt n, UInt k, Real alpha, const Real * A, const Real * B, __attribute__ ((unused)) Real beta, Real * C) { if(tr_A) { if(tr_B) matrixt_matrixt(m, n, k, A, B, C, alpha); else matrixt_matrix(m, n, k, A, B, C, alpha); } else { if(tr_B) matrix_matrixt(m, n, k, A, B, C, alpha); else matrix_matrix(m, n, k, A, B, C, alpha); } } /* -------------------------------------------------------------------------- */ template inline void Math::matVectMul(UInt m, UInt n, Real alpha, const Real * A, const Real * x, __attribute__ ((unused)) Real beta, Real * y) { if(tr_A) { matrixt_vector(m, n, A, x, y, alpha); } else { matrix_vector(m, n, A, x, y, alpha); } } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_LAPACK inline void Math::matrixEig(UInt n, Real * A, Real * d, Real * V) { // Matrix A is row major, so the lapack function in fortran will process // A^t. Asking for the left eigenvectors of A^t will give the transposed right // eigenvectors of A so in the C++ code the right eigenvectors. char jobvl; if(V != NULL) jobvl = 'V'; // compute left eigenvectors else jobvl = 'N'; // compute left eigenvectors char jobvr('N'); // compute right eigenvectors double * di = new double[n]; // imaginary part of the eigenvalues int info; int N = n; double wkopt; int lwork = -1; // query and allocate the optimal workspace dgeev_(&jobvl, &jobvr, &N, A, &N, d, di, V, &N, NULL, &N, &wkopt, &lwork, &info); lwork = int(wkopt); double * work = new double[lwork]; // solve the eigenproblem dgeev_(&jobvl, &jobvr, &N, A, &N, d, di, V, &N, NULL, &N, work, &lwork, &info); AKANTU_DEBUG_ASSERT(info == 0, "Problem computing eigenvalues/vectors. DGEEV exited with the value " << info); delete [] work; delete [] di; // I hope for you that there was no complex eigenvalues !!! } #else inline void Math::matrixEig(__attribute__((unused)) UInt n, __attribute__((unused)) Real * A, __attribute__((unused)) Real * d, __attribute__((unused)) Real * V) { AKANTU_DEBUG_ERROR("You have to compile with the support of LAPACK activated to use this function!"); } #endif /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_LAPACK inline void Math::inv(UInt n, const Real * A, Real * invA) { int N = n; int info; int * ipiv = new int[N+1]; int lwork = N*N; double * work = new double[lwork]; std::copy(A, A + n*n, invA); dgetrf_(&N, &N, invA, &N, ipiv, &info); if(info > 0) { AKANTU_DEBUG_ERROR("Singular matrix - cannot factorize it (info: " << info <<" )"); } dgetri_(&N, invA, &N, ipiv, work, &lwork, &info); if(info != 0) { AKANTU_DEBUG_ERROR("Cannot invert the matrix (info: "<< info <<" )"); } delete [] ipiv; delete [] work; } #else inline void Math::inv(__attribute__((unused)) UInt n, __attribute__((unused)) const Real * A, __attribute__((unused)) Real * Ainv) { AKANTU_DEBUG_ERROR("You have to compile with the support of LAPACK activated to use this function!"); } #endif /* -------------------------------------------------------------------------- */ inline void Math::matrix22_eigenvalues(Real * A, Real *Adiag) { ///d = determinant of Matrix A Real d = det2(A); ///b = trace of Matrix A Real b = A[0]+A[3]; Real c = sqrt(b*b - 4 *d); Adiag[0]= .5*(b + c); Adiag[1]= .5*(b - c); } /* -------------------------------------------------------------------------- */ inline void Math::matrix33_eigenvalues(Real * A, Real *Adiag) { /// a L^3 + b L^2 + c L + d = 0 matrixEig(3, A, Adiag); // Real a = -1 ; // ///b = trace of Matrix A // Real b = A[0]+A[4]+A[8]; // /// c = 0.5*(trace(M^2)-trace(M)^2) // Real c = A[1]*A[3] + A[2]*A[6] + A[5]*A[7] - A[0]*A[4] - // A[0]*A[8] - A[4]*A[8]; // ///d = determinant of Matrix A // Real d = det3(A); // // /// Define x, y, z // Real x = c/a - b*b/(3.*a*a); // Real y = 2.*b*b*b/(27.*a*a*a) - b*c/(3.*a*a) + d/a; // Real z = y*y/4. + x*x*x/27.; // /// Define I, j, k, m, n, p (so equations are not so cluttered) // Real i = sqrt(y*y/4. - z); // Real j = pow(i,1./3.); // Real k = 0; // if (std::abs(i) > 1e-12) // k = acos(-(y/(2.*i))); // // Real m = cos(k/3); // Real n = sqrt(3.)*sin(k/3); // Real p = -b/(3.*a); // // Adiag[0] = 2*j*m + p; // Adiag[1] = -j *(m + n) + p; // Adiag[2] = -j * (m - n) + p; } /* -------------------------------------------------------------------------- */ template inline void Math::eigenvalues(Real * A, Real * d) { if(dim == 1) { d[0] = A[0]; } else if(dim == 2) { matrix22_eigenvalues(A, d); } // else if(dim == 3) { matrix33_eigenvalues(A, d); } else matrixEig(dim, A, d); } /* -------------------------------------------------------------------------- */ inline Real Math::det2(const Real * mat) { return mat[0]*mat[3] - mat[1]*mat[2]; } /* -------------------------------------------------------------------------- */ inline Real Math::det3(const Real * mat) { return mat[0]*(mat[4]*mat[8]-mat[7]*mat[5]) - mat[3]*(mat[1]*mat[8]-mat[7]*mat[2]) + mat[6]*(mat[1]*mat[5]-mat[4]*mat[2]); } /* -------------------------------------------------------------------------- */ inline void Math::normal2(const Real * vec,Real * normal) { normal[0] = vec[1]; normal[1] = -vec[0]; Math::normalize2(normal); } /* -------------------------------------------------------------------------- */ inline void Math::normal3(const Real * vec1,const Real * vec2,Real * normal) { Math::vectorProduct3(vec1,vec2,normal); Math::normalize3(normal); } /* -------------------------------------------------------------------------- */ inline void Math::normalize2(Real * vec) { Real norm = Math::norm2(vec); vec[0] /= norm; vec[1] /= norm; } /* -------------------------------------------------------------------------- */ inline void Math::normalize3(Real * vec) { Real norm = Math::norm3(vec); vec[0] /= norm; vec[1] /= norm; vec[2] /= norm; } /* -------------------------------------------------------------------------- */ inline Real Math::norm2(const Real * vec) { return sqrt(vec[0]*vec[0] + vec[1]*vec[1]); } /* -------------------------------------------------------------------------- */ inline Real Math::norm3(const Real * vec) { return sqrt(vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2]); } /* -------------------------------------------------------------------------- */ inline Real Math::norm(UInt n, const Real * vec) { Real norm = 0.; for (UInt i = 0; i < n; ++i) { norm += vec[i]*vec[i]; } return sqrt(norm); } /* -------------------------------------------------------------------------- */ inline void Math::inv2(const Real * mat,Real * inv) { Real det_mat = det2(mat); inv[0] = mat[3] / det_mat; inv[1] = -mat[1] / det_mat; inv[2] = -mat[2] / det_mat; inv[3] = mat[0] / det_mat; } /* -------------------------------------------------------------------------- */ inline void Math::inv3(const Real * mat,Real * inv) { Real det_mat = det3(mat); inv[0] = (mat[4]*mat[8] - mat[7]*mat[5])/det_mat; inv[1] = (mat[2]*mat[7] - mat[8]*mat[1])/det_mat; inv[2] = (mat[1]*mat[5] - mat[4]*mat[2])/det_mat; inv[3] = (mat[5]*mat[6] - mat[8]*mat[3])/det_mat; inv[4] = (mat[0]*mat[8] - mat[6]*mat[2])/det_mat; inv[5] = (mat[2]*mat[3] - mat[5]*mat[0])/det_mat; inv[6] = (mat[3]*mat[7] - mat[6]*mat[4])/det_mat; inv[7] = (mat[1]*mat[6] - mat[7]*mat[0])/det_mat; inv[8] = (mat[0]*mat[4] - mat[3]*mat[1])/det_mat; } /* -------------------------------------------------------------------------- */ inline void Math::vectorProduct3(const Real * v1, const Real * v2, Real * res) { res[0] = v1[1]*v2[2] - v1[2]*v2[1]; res[1] = v1[2]*v2[0] - v1[0]*v2[2]; res[2] = v1[0]*v2[1] - v1[1]*v2[0]; } /* -------------------------------------------------------------------------- */ inline Real Math::vectorDot2(const Real * v1, const Real * v2) { return (v1[0]*v2[0] + v1[1]*v2[1]); } /* -------------------------------------------------------------------------- */ inline Real Math::vectorDot3(const Real * v1, const Real * v2) { return (v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]); } /* -------------------------------------------------------------------------- */ inline Real Math::distance_2d(const Real * x, const Real * y) { return sqrt((y[0] - x[0])*(y[0] - x[0]) + (y[1] - x[1])*(y[1] - x[1])); } /* -------------------------------------------------------------------------- */ inline Real Math::triangle_inradius(const Real * coord1, const Real * coord2, const Real * coord3) { /** * @f{eqnarray*}{ * r &=& A / s \\ * A &=& 1/4 * \sqrt{(a + b + c) * (a - b + c) * (a + b - c) (-a + b + c)} \\ * s &=& \frac{a + b + c}{2} * @f} */ Real a, b, c; a = distance_2d(coord1, coord2); b = distance_2d(coord2, coord3); c = distance_2d(coord1, coord3); Real s; s = (a + b + c) * 0.5; return sqrt((s - a) * (s - b) * (s - c) / s); } /* -------------------------------------------------------------------------- */ inline Real Math::distance_3d(const Real * x, const Real * y) { return sqrt((y[0] - x[0])*(y[0] - x[0]) + (y[1] - x[1])*(y[1] - x[1]) + (y[2] - x[2])*(y[2] - x[2]) ); } /* -------------------------------------------------------------------------- */ inline Real Math::tetrahedron_volume(const Real * coord1, const Real * coord2, const Real * coord3, const Real * coord4) { Real xx[9], vol; xx[0] = coord2[0]; xx[1] = coord2[1]; xx[2] = coord2[2]; xx[3] = coord3[0]; xx[4] = coord3[1]; xx[5] = coord3[2]; xx[6] = coord4[0]; xx[7] = coord4[1]; xx[8] = coord4[2]; vol = det3(xx); xx[0] = coord1[0]; xx[1] = coord1[1]; xx[2] = coord1[2]; xx[3] = coord3[0]; xx[4] = coord3[1]; xx[5] = coord3[2]; xx[6] = coord4[0]; xx[7] = coord4[1]; xx[8] = coord4[2]; vol -= det3(xx); xx[0] = coord1[0]; xx[1] = coord1[1]; xx[2] = coord1[2]; xx[3] = coord2[0]; xx[4] = coord2[1]; xx[5] = coord2[2]; xx[6] = coord4[0]; xx[7] = coord4[1]; xx[8] = coord4[2]; vol += det3(xx); xx[0] = coord1[0]; xx[1] = coord1[1]; xx[2] = coord1[2]; xx[3] = coord2[0]; xx[4] = coord2[1]; xx[5] = coord2[2]; xx[6] = coord3[0]; xx[7] = coord3[1]; xx[8] = coord3[2]; vol -= det3(xx); vol /= 6; return vol; } /* -------------------------------------------------------------------------- */ inline Real Math::tetrahedron_inradius(const Real * coord1, const Real * coord2, const Real * coord3, const Real * coord4) { Real l12, l13, l14, l23, l24, l34; l12 = distance_3d(coord1, coord2); l13 = distance_3d(coord1, coord3); l14 = distance_3d(coord1, coord4); l23 = distance_3d(coord2, coord3); l24 = distance_3d(coord2, coord4); l34 = distance_3d(coord3, coord4); Real s1, s2, s3, s4; s1 = (l12 + l23 + l13) * 0.5; s1 = sqrt(s1*(s1-l12)*(s1-l23)*(s1-l13)); s2 = (l12 + l24 + l14) * 0.5; s2 = sqrt(s2*(s2-l12)*(s2-l24)*(s2-l14)); s3 = (l23 + l34 + l24) * 0.5; s3 = sqrt(s3*(s3-l23)*(s3-l34)*(s3-l24)); s4 = (l13 + l34 + l14) * 0.5; s4 = sqrt(s4*(s4-l13)*(s4-l34)*(s4-l14)); Real volume = Math::tetrahedron_volume(coord1,coord2,coord3,coord4); return 3*volume/(s1+s2+s3+s4); } /* -------------------------------------------------------------------------- */ inline void Math::barycenter(const Real * coord, UInt nb_points, UInt spatial_dimension, Real * barycenter) { memset(barycenter, 0, spatial_dimension * sizeof(Real)); for (UInt n = 0; n < nb_points; ++n) { UInt offset = n * spatial_dimension; for (UInt i = 0; i < spatial_dimension; ++i) { barycenter[i] += coord[offset + i] / (Real) nb_points; } } } /* -------------------------------------------------------------------------- */ inline void Math::vector_2d(const Real * x, const Real * y, Real * res) { res[0] = y[0]-x[0]; res[1] = y[1]-x[1]; } /* -------------------------------------------------------------------------- */ inline void Math::vector_3d(const Real * x, const Real * y, Real * res) { res[0] = y[0]-x[0]; res[1] = y[1]-x[1]; res[2] = y[2]-x[2]; } /* -------------------------------------------------------------------------- */ inline bool Math::are_float_equal(const Real x, const Real y){ return (std::abs( x - y) < tolerance); } /* -------------------------------------------------------------------------- */ inline bool Math::isnan(Real x) { #if defined(__INTEL_COMPILER) #pragma warning ( push ) #pragma warning ( disable : 1572 ) #endif //defined(__INTEL_COMPILER) // x = x return false means x = quiet_NaN return !(x == x); #if defined(__INTEL_COMPILER) #pragma warning ( pop ) #endif //defined(__INTEL_COMPILER) } /* -------------------------------------------------------------------------- */ inline bool Math::are_vector_equal(UInt n, Real * x, Real * y){ bool test = true; for (UInt i = 0; i < n; ++i) { test &= are_float_equal(x[i],y[i]); } return test; } /* -------------------------------------------------------------------------- */ inline bool Math::intersects(Real x_min, Real x_max, Real y_min, Real y_max) { return ! ((x_max < y_min) || (x_min > y_max)); } /* -------------------------------------------------------------------------- */ inline bool Math::is_in_range(Real a, Real x_min, Real x_max) { return ((a >= x_min) && (a <= x_max)); } diff --git a/src/common/aka_vector_tmpl.hh b/src/common/aka_vector_tmpl.hh index 9c88ec7f6..2820b2a56 100644 --- a/src/common/aka_vector_tmpl.hh +++ b/src/common/aka_vector_tmpl.hh @@ -1,846 +1,846 @@ /** * @file aka_vector_inline_impl.cc * @author Nicolas Richart * @date Thu Jul 15 00:09:33 2010 * * @brief Inline functions of the classes Vector and VectorBase * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ /* Inline Functions Vector */ /* -------------------------------------------------------------------------- */ __END_AKANTU__ #include __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template inline T & Vector::operator()(UInt i, UInt j) { AKANTU_DEBUG_ASSERT(size > 0, "The vector \"" << id << "\" is empty"); AKANTU_DEBUG_ASSERT((i < size) && (j < nb_component), "The value at position [" << i << "," << j << "] is out of range in vector \"" << id << "\""); return values[i*nb_component + j]; } /* -------------------------------------------------------------------------- */ template inline const T & Vector::operator()(UInt i, UInt j) const { AKANTU_DEBUG_ASSERT(size > 0, "The vector \"" << id << "\" is empty"); AKANTU_DEBUG_ASSERT((i < size) && (j < nb_component), "The value at position [" << i << "," << j << "] is out of range in vector \"" << id << "\""); return values[i*nb_component + j]; } /* -------------------------------------------------------------------------- */ template inline void Vector::push_back(const T & value) { UInt pos = size; resizeUnitialized(size+1); std::uninitialized_fill_n(values + pos * nb_component, nb_component, value); } /* -------------------------------------------------------------------------- */ template inline void Vector::push_back(const T new_elem[]) { UInt pos = size; resizeUnitialized(size+1); T * tmp = values + nb_component * pos; std::uninitialized_copy(new_elem, new_elem + nb_component, tmp); } /* -------------------------------------------------------------------------- */ template template inline void Vector::push_back(const Vector::iterator & it) { UInt pos = size; resizeUnitialized(size+1); T * tmp = values + nb_component * pos; T * new_elem = it.data(); std::uninitialized_copy(new_elem, new_elem + nb_component, tmp); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template inline void Vector::erase(UInt i){ AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT((size > 0), "The vector is empty"); AKANTU_DEBUG_ASSERT((i < size), "The element at position [" << i << "] is out of range (" << i << ">=" << size << ")"); if(i != (size - 1)) { for (UInt j = 0; j < nb_component; ++j) { values[i*nb_component + j] = values[(size-1)*nb_component + j]; } } resize(size - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template Vector & Vector::operator-=(const Vector & vect) { AKANTU_DEBUG_ASSERT((size == vect.size) && (nb_component == vect.nb_component), "The too vector don't have the same sizes"); T * a = values; T * b = vect.values; for (UInt i = 0; i < size*nb_component; ++i) { *a -= *b; ++a;++b; } return *this; } /* -------------------------------------------------------------------------- */ template Vector & Vector::operator+=(const Vector & vect) { AKANTU_DEBUG_ASSERT((size == vect.size) && (nb_component == vect.nb_component), "The too vector don't have the same sizes"); T * a = values; T * b = vect.values; for (UInt i = 0; i < size*nb_component; ++i) { *a++ += *b++; } return *this; } /* -------------------------------------------------------------------------- */ template Vector & Vector::operator*=(const T & alpha) { T * a = values; for (UInt i = 0; i < size*nb_component; ++i) { *a++ *= alpha; } return *this; } /* -------------------------------------------------------------------------- */ /* Functions Vector */ /* -------------------------------------------------------------------------- */ template Vector::Vector (UInt size, UInt nb_component, const ID & id) : VectorBase(id), values(NULL) { AKANTU_DEBUG_IN(); allocate(size, nb_component); if(!is_scal) { T val = T(); std::uninitialized_fill(values, values + size*nb_component, val); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template Vector::Vector (UInt size, UInt nb_component, const T def_values[], const ID & id) : VectorBase(id), values(NULL) { AKANTU_DEBUG_IN(); allocate(size, nb_component); T * tmp = values; for (UInt i = 0; i < size; ++i) { tmp = values + nb_component * i; std::uninitialized_copy(def_values, def_values + nb_component, tmp); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template Vector::Vector (UInt size, UInt nb_component, const T & value, const ID & id) : VectorBase(id), values(NULL) { AKANTU_DEBUG_IN(); allocate(size, nb_component); std::uninitialized_fill_n(values, size*nb_component, value); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template Vector::Vector(const Vector & vect, bool deep, const ID & id) { AKANTU_DEBUG_IN(); this->id = (id == "") ? vect.id : id; if (deep) { allocate(vect.size, vect.nb_component); T * tmp = values; std::uninitialized_copy(vect.values, vect.values + size * nb_component, tmp); } else { this->values = vect.values; this->size = vect.size; this->nb_component = vect.nb_component; this->allocated_size = vect.allocated_size; this->size_of_type = vect.size_of_type; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template Vector::Vector(const std::vector& vect) { AKANTU_DEBUG_IN(); this->id = ""; allocate(vect.size(), 1); T * tmp = values; std::uninitialized_copy(&(vect[0]), &(vect[size-1]), tmp); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template Vector::~Vector () { AKANTU_DEBUG_IN(); AKANTU_DEBUG(dblAccessory, "Freeing " << allocated_size*nb_component*sizeof(T) / 1024. << "kB (" << id <<")"); if(values){ if(!is_scal) for (UInt i = 0; i < size * nb_component; ++i) { T * obj = values+i; obj->~T(); } free(values); } size = allocated_size = 0; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Vector::allocate(UInt size, UInt nb_component) { AKANTU_DEBUG_IN(); if (size == 0){ values = NULL; } else { values = static_cast(malloc(nb_component * size * sizeof(T))); AKANTU_DEBUG_ASSERT(values != NULL, "Cannot allocate " << nb_component * size * sizeof(T) / 1024. << "kB (" << id <<")"); } if (values == NULL) { this->size = this->allocated_size = 0; } else { AKANTU_DEBUG(dblAccessory, "Allocated " << size * nb_component * sizeof(T) / 1024. << "kB (" << id <<")"); this->size = this->allocated_size = size; } this->size_of_type = sizeof(T); this->nb_component = nb_component; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Vector::resize(UInt new_size) { UInt old_size = size; T * old_values = values; if(new_size < size) { for (UInt i = new_size * nb_component; i < size * nb_component; ++i) { T * obj = old_values+i; obj->~T(); } } resizeUnitialized(new_size); T val = T(); if(size > old_size) std::uninitialized_fill(values + old_size*nb_component, values + size*nb_component, val); } /* -------------------------------------------------------------------------- */ template void Vector::resizeUnitialized(UInt new_size) { // AKANTU_DEBUG_IN(); // free some memory if(new_size <= allocated_size) { if(allocated_size - new_size > AKANTU_MIN_ALLOCATION) { AKANTU_DEBUG(dblAccessory, "Freeing " << (allocated_size - size)*nb_component*sizeof(T) / 1024. << "kB (" << id <<")"); // Normally there are no allocation problem when reducing an array T * tmp_ptr = static_cast(realloc(values, new_size * nb_component * sizeof(T))); if(new_size != 0 && tmp_ptr == NULL) { AKANTU_DEBUG_ERROR("Cannot free data (" << id << ")" << " [current allocated size : " << allocated_size << " | " << "requested size : " << new_size << "]"); } values = tmp_ptr; allocated_size = new_size; } size = new_size; // AKANTU_DEBUG_OUT(); return; } // allocate more memory UInt size_to_alloc = (new_size - allocated_size < AKANTU_MIN_ALLOCATION) ? allocated_size + AKANTU_MIN_ALLOCATION : new_size; T *tmp_ptr = static_cast(realloc(values, size_to_alloc * nb_component * sizeof(T))); AKANTU_DEBUG_ASSERT(tmp_ptr != NULL, "Cannot allocate " << size_to_alloc * nb_component * sizeof(T) / 1024. << "kB"); if (tmp_ptr == NULL) { AKANTU_DEBUG_ERROR("Cannot allocate more data (" << id << ")" << " [current allocated size : " << allocated_size << " | " << "requested size : " << new_size << "]"); } AKANTU_DEBUG(dblAccessory, "Allocating " << (size_to_alloc - allocated_size)*nb_component*sizeof(T) / 1024. << "kB"); allocated_size = size_to_alloc; size = new_size; values = tmp_ptr; // AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Vector::extendComponentsInterlaced(UInt multiplicator, UInt block_size) { AKANTU_DEBUG_IN(); if (multiplicator == 1) return; AKANTU_DEBUG_ASSERT(multiplicator > 1, "invalid multiplicator"); AKANTU_DEBUG_ASSERT(nb_component%block_size == 0, "stride must divide actual number of components"); values = static_cast(realloc(values, nb_component*multiplicator*size* sizeof(T))); UInt new_component = nb_component/block_size * multiplicator; for (UInt i = 0,k=size-1; i < size; ++i,--k) { for (UInt j = 0; j < new_component; ++j) { UInt m = new_component - j -1; UInt n = m/multiplicator; for (UInt l = 0,p=block_size-1; l < block_size; ++l,--p) { values[k*nb_component*multiplicator+m*block_size+p] = values[k*nb_component+n*block_size+p]; } } } nb_component = nb_component * multiplicator; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template Int Vector::find(const T & elem) const { AKANTU_DEBUG_IN(); UInt i = 0; for (; (i < size) && (values[i] != elem); ++i); AKANTU_DEBUG_OUT(); return (i == size) ? -1 : (Int) i; } /* -------------------------------------------------------------------------- */ template Int Vector::find(T elem[]) const { AKANTU_DEBUG_IN(); T * it = values; UInt i = 0; for (;i < size; ++i) { if(*it == elem[0]) { T * cit = it; UInt c = 0; for(; (c < nb_component) && (*cit == elem[c]); ++c, ++cit); if(c == nb_component) { AKANTU_DEBUG_OUT(); return i; } } it += nb_component; } return -1; } /* -------------------------------------------------------------------------- */ template void Vector::copy(const Vector& vect) { AKANTU_DEBUG_IN(); if(AKANTU_DEBUG_TEST(dblWarning)) if(vect.nb_component != nb_component) { AKANTU_DEBUG(dblWarning, "The two vectors do not have the same number of components"); } // this->id = vect.id; resize((vect.size * vect.nb_component) / nb_component); T * tmp = values; std::uninitialized_copy(vect.values, vect.values + size * nb_component, tmp); // memcpy(this->values, vect.values, vect.size * vect.nb_component * sizeof(T)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Vector::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); Real real_size = allocated_size * nb_component * size_of_type / 1024.0; std::streamsize prec = stream.precision(); std::ios_base::fmtflags ff = stream.flags(); stream.setf (std::ios_base::showbase); stream.precision(2); stream << space << "Vector<" << debug::demangle(typeid(T).name()) << "> [" << std::endl; stream << space << " + id : " << this->id << std::endl; stream << space << " + size : " << this->size << std::endl; stream << space << " + nb_component : " << this->nb_component << std::endl; stream << space << " + allocated size : " << this->allocated_size << std::endl; stream << space << " + memory size : " << real_size << "kB" << std::endl; if(!AKANTU_DEBUG_LEVEL_IS_TEST()) stream << space << " + address : " << std::hex << this->values << std::dec << std::endl; stream.precision(prec); stream.flags(ff); if(AKANTU_DEBUG_TEST(dblDump)) { stream << space << " + values : {"; for (UInt i = 0; i < this->size; ++i) { stream << "{"; for (UInt j = 0; j < this->nb_component; ++j) { stream << this->values[i*nb_component + j]; if(j != this->nb_component - 1) stream << ", "; } stream << "}"; if(i != this->size - 1) stream << ", "; } stream << "}" << std::endl; } stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ /* Inline Functions VectorBase */ /* -------------------------------------------------------------------------- */ inline UInt VectorBase::getMemorySize() const { return allocated_size * nb_component * size_of_type; } inline void VectorBase::empty() { size = 0; } /* -------------------------------------------------------------------------- */ /* Iterators */ /* -------------------------------------------------------------------------- */ template template class Vector::iterator_internal { public: typedef R value_type; typedef R* pointer; typedef R& reference; typedef IR internal_value_type; typedef IR* internal_pointer; typedef std::ptrdiff_t difference_type; typedef std::random_access_iterator_tag iterator_category; public: - iterator_internal() : offset(0), initial(NULL), ret(NULL) {}; + iterator_internal() : _offset(0), initial(NULL), ret(NULL) {}; - iterator_internal(pointer_type data, UInt offset) : - offset(offset), + iterator_internal(pointer_type data, UInt _offset) : + _offset(_offset), initial(data), ret(new value_type(data)) { - AKANTU_DEBUG_ASSERT(offset == ret->size(), + AKANTU_DEBUG_ASSERT(_offset == ret->size(), "The iterator_internal is not compatible with the type " << typeid(value_type).name()); } - iterator_internal(pointer warped) : offset(warped->size()), + iterator_internal(pointer warped) : _offset(warped->size()), initial(warped->storage()), ret(const_cast(warped)) { } iterator_internal(const iterator_internal & it) { if(this != &it) { - this->offset = it.offset; + this->_offset = it._offset; this->initial = it.initial; this->ret = new internal_value_type(*it.ret); } } virtual ~iterator_internal() { delete ret; }; inline iterator_internal & operator=(const iterator_internal & it) { if(this != &it) { - this->offset = it.offset; + this->_offset = it._offset; this->initial = it.initial; if(this->ret) this->ret->shallowCopy(*it.ret); else this->ret = new internal_value_type(*it.ret); } return *this; } inline reference operator*() { return *ret; }; inline pointer operator->() { return ret; }; - inline iterator_internal & operator++() { ret->values += offset; return *this; }; - inline iterator_internal & operator--() { ret->values -= offset; return *this; }; + inline iterator_internal & operator++() { ret->values += _offset; return *this; }; + inline iterator_internal & operator--() { ret->values -= _offset; return *this; }; - inline iterator_internal & operator+=(const UInt n) { ret->values += offset * n; return *this; } - inline iterator_internal & operator-=(const UInt n) { ret->values -= offset * n; return *this; } + inline iterator_internal & operator+=(const UInt n) { ret->values += _offset * n; return *this; } + inline iterator_internal & operator-=(const UInt n) { ret->values -= _offset * n; return *this; } - inline reference operator[](const UInt n) { ret->values = initial + n*offset; return *ret; } + inline reference operator[](const UInt n) { ret->values = initial + n*_offset; return *ret; } - inline bool operator==(const iterator_internal & other) const { return (*this).ret->values == other.ret->values; } - inline bool operator!=(const iterator_internal & other) const { return (*this).ret->values != other.ret->values; } - inline bool operator <(const iterator_internal & other) const { return (*this).ret->values < other.ret->values; } - inline bool operator<=(const iterator_internal & other) const { return (*this).ret->values <= other.ret->values; } - inline bool operator> (const iterator_internal & other) const { return (*this).ret->values > other.ret->values; } - inline bool operator>=(const iterator_internal & other) const { return (*this).ret->values >= other.ret->values; } + inline bool operator==(const iterator_internal & other) const { return (*this).ret->storage() == other.ret->storage(); } + inline bool operator!=(const iterator_internal & other) const { return (*this).ret->storage() != other.ret->storage(); } + inline bool operator <(const iterator_internal & other) const { return (*this).ret->storage() < other.ret->storage(); } + inline bool operator<=(const iterator_internal & other) const { return (*this).ret->storage() <= other.ret->storage(); } + inline bool operator> (const iterator_internal & other) const { return (*this).ret->storage() > other.ret->storage(); } + inline bool operator>=(const iterator_internal & other) const { return (*this).ret->storage() >= other.ret->storage(); } inline iterator_internal operator+(difference_type n) { iterator_internal tmp(*this); tmp += n; return tmp; } inline iterator_internal operator-(difference_type n) { iterator_internal tmp(*this); tmp -= n; return tmp; } inline difference_type operator-(const iterator_internal & b) { return ret->values - b.ret->values; } - inline pointer_type data() const { - return ret->storage(); - } + inline pointer_type data() const { return ret->storage(); } + inline difference_type offset() const { return _offset; } protected: - UInt offset; + UInt _offset; pointer_type initial; internal_pointer ret; }; /* -------------------------------------------------------------------------- */ template inline Vector::iterator< types::Vector > Vector::begin(UInt n) { AKANTU_DEBUG_ASSERT(nb_component == n, "The iterator is not compatible with the type Vector(" << n<< ")"); return iterator< types::Vector >(new types::Vector(values, n)); } /* -------------------------------------------------------------------------- */ template inline Vector::iterator< types::Vector > Vector::end(UInt n) { AKANTU_DEBUG_ASSERT(nb_component == n, "The iterator is not compatible with the type Vector(" << n<< ")"); return iterator< types::Vector >(new types::Vector(values + nb_component * size, n)); } /* -------------------------------------------------------------------------- */ template inline Vector::const_iterator< types::Vector > Vector::begin(UInt n) const { AKANTU_DEBUG_ASSERT(nb_component == n, "The iterator is not compatible with the type Vector(" << n<< ")"); return const_iterator< types::Vector >(new types::Vector(values, n)); } /* -------------------------------------------------------------------------- */ template inline Vector::const_iterator< types::Vector > Vector::end(UInt n) const { AKANTU_DEBUG_ASSERT(nb_component == n, "The iterator is not compatible with the type Vector(" << n<< ")"); return const_iterator< types::Vector >(new types::Vector(values + nb_component * size, n)); } /* -------------------------------------------------------------------------- */ template inline Vector::iterator< types::Matrix > Vector::begin(UInt m, UInt n) { AKANTU_DEBUG_ASSERT(nb_component == n*m, "The iterator is not compatible with the type Matrix(" << m << "," << n<< ")"); return iterator< types::Matrix >(new types::Matrix(values, m, n)); } /* -------------------------------------------------------------------------- */ template inline Vector::iterator< types::Matrix > Vector::end(UInt m, UInt n) { AKANTU_DEBUG_ASSERT(nb_component == n*m, "The iterator is not compatible with the type Matrix(" << m << "," << n<< ")"); return iterator< types::Matrix >(new types::Matrix(values + nb_component * size, m, n)); } /* -------------------------------------------------------------------------- */ template inline Vector::const_iterator< types::Matrix > Vector::begin(UInt m, UInt n) const { AKANTU_DEBUG_ASSERT(nb_component == n*m, "The iterator is not compatible with the type Matrix(" << m << "," << n<< ")"); return const_iterator< types::Matrix >(new types::Matrix(values, m, n)); } /* -------------------------------------------------------------------------- */ template inline Vector::const_iterator< types::Matrix > Vector::end(UInt m, UInt n) const { AKANTU_DEBUG_ASSERT(nb_component == n*m, "The iterator is not compatible with the type Matrix(" << m << "," << n<< ")"); return const_iterator< types::Matrix >(new types::Matrix(values + nb_component * size, m, n)); } /* -------------------------------------------------------------------------- */ template inline Vector::iterator< types::Matrix > Vector::begin_reinterpret(UInt m, UInt n, __attribute__((unused)) UInt size, __attribute__((unused)) UInt nb_component) { AKANTU_DEBUG_ASSERT(nb_component * size == this->nb_component * this->size, "The new values for size (" << size << ") and nb_component (" << nb_component << ") are not compatible with the one of this vector"); AKANTU_DEBUG_ASSERT(nb_component == n*m, "The iterator is not compatible with the type Matrix(" << m << "," << n<< ")"); return iterator< types::Matrix >(new types::Matrix(values + nb_component * 0, m, n)); } /* -------------------------------------------------------------------------- */ template inline Vector::iterator< types::Matrix > Vector::end_reinterpret(UInt m, UInt n, UInt size, UInt nb_component) { AKANTU_DEBUG_ASSERT(nb_component * size == this->nb_component * this->size, "The new values for size (" << size << ") and nb_component (" << nb_component << ") are not compatible with the one of this vector"); AKANTU_DEBUG_ASSERT(nb_component == n*m, "The iterator is not compatible with the type Matrix(" << m << "," << n<< ")"); return iterator< types::Matrix >(new types::Matrix(values + nb_component * size, m, n)); } /* -------------------------------------------------------------------------- */ template inline Vector::const_iterator< types::Matrix > Vector::begin_reinterpret(UInt m, UInt n, __attribute__((unused)) UInt size, __attribute__((unused)) UInt nb_component) const { AKANTU_DEBUG_ASSERT(nb_component * size == this->nb_component * this->size, "The new values for size (" << size << ") and nb_component (" << nb_component << ") are not compatible with the one of this vector"); AKANTU_DEBUG_ASSERT(nb_component == n*m, "The iterator is not compatible with the type Matrix(" << m << "," << n<< ")"); return const_iterator< types::Matrix >(new types::Matrix(values + nb_component * 0, m, n)); } /* -------------------------------------------------------------------------- */ template inline Vector::const_iterator< types::Matrix > Vector::end_reinterpret(UInt m, UInt n, UInt size, UInt nb_component) const { AKANTU_DEBUG_ASSERT(nb_component * size == this->nb_component * this->size, "The new values for size (" << size << ") and nb_component (" << nb_component << ") are not compatible with the one of this vector"); AKANTU_DEBUG_ASSERT(nb_component == n*m, "The iterator is not compatible with the type Matrix(" << m << "," << n<< ")"); return const_iterator< types::Matrix >(new types::Matrix(values + nb_component * size, m, n)); } /* -------------------------------------------------------------------------- */ /** * Specialization for scalar types */ template template class Vector::iterator_internal { public: typedef R value_type; typedef R* pointer; typedef R& reference; typedef IR internal_value_type; typedef IR* internal_pointer; typedef std::ptrdiff_t difference_type; typedef std::random_access_iterator_tag iterator_category; public: - iterator_internal(pointer data = NULL, __attribute__ ((unused)) UInt offset = 1) : ret(data), initial(data) { }; + iterator_internal(pointer data = NULL, __attribute__ ((unused)) UInt _offset = 1) : _offset(_offset), ret(data), initial(data) { }; iterator_internal(const iterator_internal & it) { if(this != &it) { this->ret = it.ret; this->initial = it.initial; } } virtual ~iterator_internal() { }; inline iterator_internal & operator=(const iterator_internal & it) { if(this != &it) { this->ret = it.ret; this->initial = it.initial; } return *this; } inline reference operator*() { return *ret; }; inline pointer operator->() { return ret; }; inline iterator_internal & operator++() { ++ret; return *this; }; inline iterator_internal & operator--() { --ret; return *this; }; inline iterator_internal & operator+=(const UInt n) { ret += n; return *this; } inline iterator_internal & operator-=(const UInt n) { ret -= n; return *this; } inline reference operator[](const UInt n) { ret = initial + n; return *ret; } inline bool operator==(const iterator_internal & other) const { return ret == other.ret; } inline bool operator!=(const iterator_internal & other) const { return ret != other.ret; } inline bool operator< (const iterator_internal & other) const { return ret < other.ret; } inline bool operator<=(const iterator_internal & other) const { return ret <= other.ret; } inline bool operator> (const iterator_internal & other) const { return ret > other.ret; } inline bool operator>=(const iterator_internal & other) const { return ret >= other.ret; } inline iterator_internal operator-(difference_type n) { return iterator_internal(ret - n); } inline iterator_internal operator+(difference_type n) { return iterator_internal(ret + n); } inline difference_type operator-(const iterator_internal & b) { return ret - b.ret; } inline pointer data() const { return ret; } - + inline difference_type offset() const { return _offset; } protected: + difference_type _offset; pointer ret; pointer initial; }; /* -------------------------------------------------------------------------- */ template template inline void Vector::erase(const iterator & it) { T * curr = it.getCurrentStorage(); UInt pos = (curr - values) / nb_component; erase(pos); } /* -------------------------------------------------------------------------- */ template inline Vector::iterator Vector::begin() { AKANTU_DEBUG_ASSERT(nb_component == 1, "this iterator cannot be used on a vector which has nb_component != 1"); return iterator(values); } /* -------------------------------------------------------------------------- */ template inline Vector::iterator Vector::end() { AKANTU_DEBUG_ASSERT(nb_component == 1, "this iterator cannot be used on a vector which has nb_component != 1"); return iterator(values + size); } /* -------------------------------------------------------------------------- */ template inline Vector::const_iterator Vector::begin() const { AKANTU_DEBUG_ASSERT(nb_component == 1, "this iterator cannot be used on a vector which has nb_component != 1"); return const_iterator(values); } /* -------------------------------------------------------------------------- */ template inline Vector::const_iterator Vector::end() const { AKANTU_DEBUG_ASSERT(nb_component == 1, "this iterator cannot be used on a vector which has nb_component != 1"); return const_iterator(values + size); } diff --git a/src/fem/by_element_type.hh b/src/fem/by_element_type.hh index 430535cc5..7afb61d73 100644 --- a/src/fem/by_element_type.hh +++ b/src/fem/by_element_type.hh @@ -1,185 +1,186 @@ /** * @file by_element_type.hh * @author Nicolas Richart * @date Thu Aug 4 14:40:34 2011 * * @brief storage class by element type * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_BY_ELEMENT_TYPE_HH__ #define __AKANTU_BY_ELEMENT_TYPE_HH__ #include "aka_common.hh" #include "aka_vector.hh" #include "aka_memory.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* ByElementType */ /* -------------------------------------------------------------------------- */ template class ByElementType { protected: typedef std::map DataMap; public: ByElementType(const ID & id = "by_element_type", const ID & parent_id = ""); ~ByElementType(); inline static std::string printType(const ElementType & type, const GhostType & ghost_type); inline bool exists(ElementType type, GhostType ghost_type = _not_ghost) const; inline const Stored & operator()(const ElementType & type, const GhostType & ghost_type = _not_ghost) const; inline Stored & operator()(const ElementType & type, const GhostType & ghost_type = _not_ghost); inline Stored & operator()(const Stored & insert, const ElementType & type, const GhostType & ghost_type = _not_ghost); void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* Element type Iterator */ /* ------------------------------------------------------------------------ */ class type_iterator : private std::iterator { public: typedef const ElementType value_type; typedef const ElementType* pointer; typedef const ElementType& reference; protected: typedef typename ByElementType::DataMap::const_iterator DataMapIterator; public: type_iterator(DataMapIterator & list_begin, DataMapIterator & list_end, UInt dim, ElementKind ek); type_iterator(const type_iterator & it); inline reference operator*(); + inline const reference operator*() const; inline type_iterator & operator++(); type_iterator operator++(int); - inline bool operator==(const type_iterator & other); - inline bool operator!=(const type_iterator & other); + inline bool operator==(const type_iterator & other) const; + inline bool operator!=(const type_iterator & other) const; private: DataMapIterator list_begin; DataMapIterator list_end; UInt dim; ElementKind kind; }; inline type_iterator firstType(UInt dim = 0, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const; inline type_iterator lastType(UInt dim = 0, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const; inline void setID(const ID & id) { this->id = id; } protected: inline DataMap & getData(GhostType ghost_type); inline const DataMap & getData(GhostType ghost_type) const; public: AKANTU_GET_MACRO(ID, id, ID); /* -------------------------------------------------------------------------- */ protected: ID id; DataMap data; DataMap ghost_data; }; /* -------------------------------------------------------------------------- */ /* Some typedefs */ /* -------------------------------------------------------------------------- */ template class ByElementTypeVector : public ByElementType *>, protected Memory { protected: typedef typename ByElementType *>::DataMap DataMap; public: ByElementTypeVector() {}; // ByElementTypeVector(const ID & id = "by_element_type_vector", // const MemoryID & memory_id = 0) : // ByElementType *>(id, memory_id) {}; ByElementTypeVector(const ID & id, const ID & parent_id, const MemoryID & memory_id = 0) : ByElementType *>(id, parent_id), Memory(memory_id) {}; inline Vector & alloc(UInt size, UInt nb_component, const ElementType & type, const GhostType & ghost_type); inline void alloc(UInt size, UInt nb_component, const ElementType & type); inline const Vector & operator()(const ElementType & type, const GhostType & ghost_type = _not_ghost) const; inline Vector & operator()(const ElementType & type, const GhostType & ghost_type = _not_ghost); inline void setVector(const ElementType & type, const GhostType & ghost_type, const Vector & vect); inline void free(); inline void onElementsRemoved(const ByElementTypeVector & new_numbering); }; /// to store data Vector by element type typedef ByElementTypeVector ByElementTypeReal; /// to store data Vector by element type typedef ByElementTypeVector ByElementTypeInt; /// to store data Vector by element type typedef ByElementTypeVector ByElementTypeUInt; /// Map of data of type UInt stored in a mesh typedef std::map *> UIntDataMap; typedef ByElementType ByElementTypeUIntDataMap; // /* -------------------------------------------------------------------------- */ // /* inline functions */ // /* -------------------------------------------------------------------------- */ // #if defined (AKANTU_INCLUDE_INLINE_IMPL) // # include "by_element_type_inline_impl.cc" // #endif __END_AKANTU__ #endif /* __AKANTU_BY_ELEMENT_TYPE_HH__ */ diff --git a/src/fem/by_element_type_tmpl.hh b/src/fem/by_element_type_tmpl.hh index 6423556f0..e2356f734 100644 --- a/src/fem/by_element_type_tmpl.hh +++ b/src/fem/by_element_type_tmpl.hh @@ -1,399 +1,406 @@ /** * @file by_element_type_inline_impl.cc * @author Nicolas Richart * @date Thu Aug 4 14:41:29 2011 * * @brief * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* ByElementType */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template inline std::string ByElementType::printType(const ElementType & type, const GhostType & ghost_type) { std::stringstream sstr; sstr << "(" << ghost_type << ":" << type << ")"; return sstr.str(); } /* -------------------------------------------------------------------------- */ template inline bool ByElementType::exists(ElementType type, GhostType ghost_type) const { return this->getData(ghost_type).find(type) != this->getData(ghost_type).end(); } /* -------------------------------------------------------------------------- */ template inline const Stored & ByElementType::operator()(const ElementType & type, const GhostType & ghost_type) const { typename DataMap::const_iterator it = this->getData(ghost_type).find(type); if(it == this->getData(ghost_type).end()) AKANTU_EXCEPTION("No element of type " << ByElementType::printType(type, ghost_type) << " in this ByElementType<" << debug::demangle(typeid(Stored).name()) << "> class"); return it->second; } /* -------------------------------------------------------------------------- */ template inline Stored & ByElementType::operator()(const ElementType & type, const GhostType & ghost_type) { typename DataMap::iterator it = this->getData(ghost_type).find(type); // if(it == this->getData(ghost_type).end()) // AKANTU_EXCEPTION("No element of type " // << ByElementType::printType(type, ghost_type) // << " in this ByElementType<" // << debug::demangle(typeid(Stored).name()) << "> class"); if(it == this->getData(ghost_type).end()) { DataMap & data = this->getData(ghost_type); const std::pair & res = data.insert(std::pair(type, Stored())); it = res.first; } return it->second; } /* -------------------------------------------------------------------------- */ template inline Stored & ByElementType::operator()(const Stored & insert, const ElementType & type, const GhostType & ghost_type) { typename DataMap::iterator it = this->getData(ghost_type).find(type); if(it != this->getData(ghost_type).end()) { AKANTU_EXCEPTION("Element of type " << ByElementType::printType(type, ghost_type) << " already in this ByElementType<" << debug::demangle(typeid(Stored).name()) << "> class"); } else { DataMap & data = this->getData(ghost_type); const std::pair & res = data.insert(std::pair(type, insert)); it = res.first; } return it->second; } /* -------------------------------------------------------------------------- */ template inline typename ByElementType::DataMap & ByElementType::getData(GhostType ghost_type) { if(ghost_type == _not_ghost) return data; else return ghost_data; } /* -------------------------------------------------------------------------- */ template inline const typename ByElementType::DataMap & ByElementType::getData(GhostType ghost_type) const { if(ghost_type == _not_ghost) return data; else return ghost_data; } /* -------------------------------------------------------------------------- */ /// Works only if stored is a pointer to a class with a printself method template void ByElementType::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "ByElementType<" << debug::demangle(typeid(Stored).name()) << "> [" << std::endl; for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; const DataMap & data = getData(gt); typename DataMap::const_iterator it; for(it = data.begin(); it != data.end(); ++it) { stream << space << space << ByElementType::printType(it->first, gt) << " [" << std::endl; it->second->printself(stream, indent + 3); } } stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ template ByElementType::ByElementType(const ID & id, const ID & parent_id) { AKANTU_DEBUG_IN(); std::stringstream sstr; if(parent_id != "") sstr << parent_id << ":"; sstr << id; this->id = sstr.str(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template ByElementType::~ByElementType() { } /* -------------------------------------------------------------------------- */ template inline Vector & ByElementTypeVector::alloc(UInt size, UInt nb_component, const ElementType & type, const GhostType & ghost_type) { std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; Vector * tmp; typename ByElementTypeVector::DataMap::iterator it = this->getData(ghost_type).find(type); if(it == this->getData(ghost_type).end()) { std::stringstream sstr; sstr << this->id << ":" << type << ghost_id; tmp = &(Memory::alloc(sstr.str(), size, nb_component, T())); std::stringstream sstrg; sstrg << ghost_type; tmp->setTag(sstrg.str()); this->getData(ghost_type)[type] = tmp; } else { AKANTU_DEBUG_INFO("The vector " << this->id << this->printType(type, ghost_type) << " already exists, it is resized instead of allocated."); tmp = it->second; it->second->resize(size); } return *tmp; } /* -------------------------------------------------------------------------- */ template inline void ByElementTypeVector::alloc(UInt size, UInt nb_component, const ElementType & type) { this->alloc(size, nb_component, type, _not_ghost); this->alloc(size, nb_component, type, _ghost); } /* -------------------------------------------------------------------------- */ template inline void ByElementTypeVector::free() { AKANTU_DEBUG_IN(); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; const DataMap & data = this->getData(gt); typename DataMap::const_iterator it; for(it = data.begin(); it != data.end(); ++it) { dealloc(it->second->getID()); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline const Vector & ByElementTypeVector::operator()(const ElementType & type, const GhostType & ghost_type) const { typename ByElementTypeVector::DataMap::const_iterator it = this->getData(ghost_type).find(type); if(it == this->getData(ghost_type).end()) AKANTU_EXCEPTION("No element of type " << ByElementTypeVector::printType(type, ghost_type) << " in this const ByElementTypeVector<" << debug::demangle(typeid(T).name()) << "> class(\"" << this->id << "\")"); return *(it->second); } /* -------------------------------------------------------------------------- */ template inline Vector & ByElementTypeVector::operator()(const ElementType & type, const GhostType & ghost_type) { typename ByElementTypeVector::DataMap::iterator it = this->getData(ghost_type).find(type); if(it == this->getData(ghost_type).end()) AKANTU_EXCEPTION("No element of type " << ByElementTypeVector::printType(type, ghost_type) << " in this ByElementTypeVector<" << debug::demangle(typeid(T).name()) << "> class (\"" << this->id << "\")"); return *(it->second); } /* -------------------------------------------------------------------------- */ template inline void ByElementTypeVector::setVector(const ElementType & type, const GhostType & ghost_type, const Vector & vect) { typename ByElementTypeVector::DataMap::iterator it = this->getData(ghost_type).find(type); if(AKANTU_DEBUG_TEST(dblWarning) && it != this->getData(ghost_type).end() && it->second != &vect) { AKANTU_DEBUG_WARNING("The Vector " << this->printType(type, ghost_type) << " is already registred, this call can lead to a memory leak."); } this->getData(ghost_type)[type] = &(const_cast &>(vect)); } /* -------------------------------------------------------------------------- */ template inline void ByElementTypeVector::onElementsRemoved(const ByElementTypeVector & new_numbering) { for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; ByElementTypeVector::type_iterator it = new_numbering.firstType(0, gt, _ek_not_defined); ByElementTypeVector::type_iterator end = new_numbering.lastType(0, gt, _ek_not_defined); for (; it != end; ++it) { ElementType type = *it; if(this->exists(type, gt)){ const Vector & renumbering = new_numbering(type, gt); Vector & vect = this->operator()(type, gt); UInt nb_component = vect.getNbComponent(); Vector tmp(renumbering.getSize(), nb_component); UInt new_size = 0; for (UInt i = 0; i < vect.getSize(); ++i) { UInt new_i = renumbering(i); if(new_i != UInt(-1)) { memcpy(tmp.storage() + new_i * nb_component, vect.storage() + i *nb_component, nb_component * sizeof(T)); ++new_size; } } tmp.resize(new_size); vect.copy(tmp); } } } } /* -------------------------------------------------------------------------- */ /* ElementType Iterator */ /* -------------------------------------------------------------------------- */ template ByElementType::type_iterator::type_iterator(DataMapIterator & list_begin, DataMapIterator & list_end, UInt dim, ElementKind ek) : list_begin(list_begin), list_end(list_end), dim(dim), kind(ek) { } /* -------------------------------------------------------------------------- */ template ByElementType::type_iterator::type_iterator(const type_iterator & it) : list_begin(it.list_begin), list_end(it.list_end), dim(it.dim), kind(it.kind) { } /* -------------------------------------------------------------------------- */ template inline typename ByElementType::type_iterator::reference ByElementType::type_iterator::operator*() { return list_begin->first; } +/* -------------------------------------------------------------------------- */ +template +inline typename ByElementType::type_iterator::reference +ByElementType::type_iterator::operator*() const { + return list_begin->first; +} + /* -------------------------------------------------------------------------- */ template inline typename ByElementType::type_iterator & ByElementType::type_iterator::operator++() { ++list_begin; while((list_begin != list_end) && (((dim != 0) && (dim != Mesh::getSpatialDimension(list_begin->first))) || ((kind != _ek_not_defined) && (kind != Mesh::getKind(list_begin->first))))) ++list_begin; return *this; } /* -------------------------------------------------------------------------- */ template typename ByElementType::type_iterator ByElementType::type_iterator::operator++(int) { type_iterator tmp(*this); operator++(); return tmp; } /* -------------------------------------------------------------------------- */ template -inline bool ByElementType::type_iterator::operator==(const type_iterator & other) { +inline bool ByElementType::type_iterator::operator==(const type_iterator & other) const { return this->list_begin == other.list_begin; } /* -------------------------------------------------------------------------- */ template -inline bool ByElementType::type_iterator::operator!=(const type_iterator & other) { +inline bool ByElementType::type_iterator::operator!=(const type_iterator & other) const { return this->list_begin != other.list_begin; } /* -------------------------------------------------------------------------- */ template inline typename ByElementType::type_iterator ByElementType::firstType(UInt dim, GhostType ghost_type, ElementKind kind) const { typename DataMap::const_iterator b,e; if(ghost_type == _not_ghost) { b = data.begin(); e = data.end(); } else { b = ghost_data.begin(); e = ghost_data.end(); } // loop until the first valid type while((b != e) && (((dim != 0) && (dim != Mesh::getSpatialDimension(b->first))) || ((kind != _ek_not_defined) && (kind != Mesh::getKind(b->first))))) ++b; return typename ByElementType::type_iterator(b, e, dim, kind); } /* -------------------------------------------------------------------------- */ template inline typename ByElementType::type_iterator ByElementType::lastType(UInt dim, GhostType ghost_type, ElementKind kind) const { typename DataMap::const_iterator e; if(ghost_type == _not_ghost) { e = data.end(); } else { e = ghost_data.end(); } return typename ByElementType::type_iterator(e, e, dim, kind); } diff --git a/src/fem/element_class_inline_impl.cc b/src/fem/element_class_inline_impl.cc index 4fa0a1f95..a55906a25 100644 --- a/src/fem/element_class_inline_impl.cc +++ b/src/fem/element_class_inline_impl.cc @@ -1,404 +1,404 @@ /** * @file element_class_inline_impl.cc * @author Nicolas Richart * @author Guillaume Anciaux * @date Thu Jul 15 10:28:28 2010 * * @brief Implementation of the inline functions of the class element_class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ template inline Real * ElementClass::getQuadraturePoints() { return quad; } /* -------------------------------------------------------------------------- */ template inline UInt ElementClass::getShapeSize() { return nb_nodes_per_element; } /* -------------------------------------------------------------------------- */ template inline UInt ElementClass::getShapeDerivativesSize() { return nb_nodes_per_element * spatial_dimension; } /* -------------------------------------------------------------------------- */ template void ElementClass::preComputeStandards(const Real * coord, const UInt dimension, Real * shape, Real * dshape, Real * jacobians) { // ask for computation of shapes computeShapes(quad, nb_quadrature_points, shape); // compute dnds Real dnds[nb_nodes_per_element * spatial_dimension * nb_quadrature_points]; computeDNDS(quad, nb_quadrature_points, dnds); // compute dxds Real dxds[dimension * spatial_dimension * nb_quadrature_points]; computeDXDS(dnds, nb_quadrature_points, coord, dimension, dxds); // jacobian computeJacobian(dxds, nb_quadrature_points, dimension, jacobians); // if dimension == spatial_dimension compute shape derivatives if (dimension == spatial_dimension) { computeShapeDerivatives(dxds, dnds, nb_quadrature_points, dimension, dshape); } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapes(const Real * natural_coords, const UInt nb_points, Real * shapes) { Real * cpoint = const_cast(natural_coords); for (UInt p = 0; p < nb_points; ++p) { computeShapes(cpoint, shapes); shapes += nb_nodes_per_element; cpoint += spatial_dimension; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapes(const Real * natural_coords, const UInt nb_points, Real * shapes, const Real * local_coord, UInt id) { Real * cpoint = const_cast(natural_coords); for (UInt p = 0; p < nb_points; ++p) { computeShapes(cpoint, shapes, local_coord, id); shapes += nb_nodes_per_element; cpoint += spatial_dimension; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeDNDS(const Real * natural_coords, const UInt nb_points, Real * dnds) { Real * cpoint = const_cast(natural_coords); Real * cdnds = dnds; for (UInt p = 0; p < nb_points; ++p) { computeDNDS(cpoint, cdnds); cpoint += spatial_dimension; cdnds += nb_nodes_per_element * spatial_dimension; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeDXDS(const Real * dnds, const UInt nb_points, const Real * node_coords, const UInt dimension, Real * dxds) { Real * cdnds = const_cast(dnds); Real * cdxds = dxds; for (UInt p = 0; p < nb_points; ++p) { computeDXDS(cdnds, node_coords, dimension, cdxds); cdnds += nb_nodes_per_element * spatial_dimension; cdxds += spatial_dimension * dimension; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeDXDS(const Real * dnds, const Real * node_coords, const UInt dimension, Real * dxds) { /// @f$ J = dxds = dnds * x @f$ Math::matrix_matrix(spatial_dimension, dimension, nb_nodes_per_element, dnds, node_coords, dxds); } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeJacobian(const Real * dxds, const UInt nb_points, const UInt dimension, Real * jac) { Real * cdxds = const_cast(dxds); Real * cjac = jac; for (UInt p = 0; p < nb_points; ++p) { computeJacobian(cdxds, dimension, *cjac); // AKANTU_DEBUG_ASSERT((cjac[0] > 0), // "Negative jacobian computed, possible problem in the element node order."); cdxds += spatial_dimension * dimension; cjac++; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapeDerivatives(const Real * dxds, const Real * dnds, const UInt nb_points, __attribute__ ((unused)) const UInt dimension, Real * shape_deriv) { AKANTU_DEBUG_ASSERT(dimension == spatial_dimension,"gradient in space " << dimension << " cannot be evaluated for element of dimension " << spatial_dimension); Real * cdxds = const_cast(dxds); Real * cdnds = const_cast(dnds); for (UInt p = 0; p < nb_points; ++p) { computeShapeDerivatives(cdxds, cdnds, shape_deriv); cdnds += spatial_dimension * nb_nodes_per_element; cdxds += spatial_dimension * spatial_dimension; shape_deriv += nb_nodes_per_element * spatial_dimension; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapeDerivatives(const Real * natural_coords, const UInt nb_points, const UInt dimension, Real * shape_deriv, const Real * local_coord, UInt id) { AKANTU_DEBUG_ASSERT(dimension == spatial_dimension,"Gradient in space " << dimension << " cannot be evaluated for element of dimension " << spatial_dimension); Real * cpoint = const_cast(natural_coords); for (UInt p = 0; p < nb_points; ++p) { computeShapeDerivatives(cpoint, shape_deriv, local_coord, id); shape_deriv += nb_nodes_per_element * dimension; cpoint += dimension; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapeDerivatives(const Real * dxds, const Real * dnds, Real * shape_deriv) { /// @f$ dxds = J^{-1} @f$ Real inv_dxds[spatial_dimension * spatial_dimension]; if (spatial_dimension == 1) inv_dxds[0] = 1./dxds[0]; if (spatial_dimension == 2) Math::inv2(dxds, inv_dxds); if (spatial_dimension == 3) Math::inv3(dxds, inv_dxds); Math::matrixt_matrixt(nb_nodes_per_element, spatial_dimension, spatial_dimension, dnds, inv_dxds, shape_deriv); } /* -------------------------------------------------------------------------- */ template inline Real ElementClass::getInradius(__attribute__ ((unused)) const Real * coord) { AKANTU_DEBUG_ERROR("Function not implemented for type : " << type); return 0; } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeNormalsOnQuadPoint(const Real * coord, const UInt dimension, Real * normals) { AKANTU_DEBUG_ASSERT((dimension - 1) == spatial_dimension, "cannot extract a normal because of dimension mismatch " << dimension << " " << spatial_dimension); Real * cpoint = const_cast(quad); Real * cnormals = normals; Real dnds[spatial_dimension*nb_nodes_per_element]; Real dxds[spatial_dimension*dimension]; for (UInt p = 0; p < nb_quadrature_points; ++p) { computeDNDS(cpoint,dnds); computeDXDS(dnds,coord,dimension,dxds); if (dimension == 2) { Math::normal2(dxds,cnormals); } if (dimension == 3){ Math::normal3(dxds,dxds+dimension,cnormals); } cpoint += spatial_dimension; cnormals += dimension; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::interpolateOnNaturalCoordinates(const Real * natural_coords, const Real * nodal_values, UInt dimension, Real * interpolated){ Real shapes[nb_nodes_per_element]; computeShapes(natural_coords,shapes); Math::matrix_matrix(1, dimension, nb_nodes_per_element, shapes, nodal_values, interpolated); } /* -------------------------------------------------------------------------- */ template inline void ElementClass::inverseMap(const types::RVector & real_coords, const types::RMatrix & node_coords, UInt dimension, types::RVector & natural_coords, Real tolerance){ //matric copy of the real_coords types::RMatrix mreal_coords(real_coords.storage(),1,spatial_dimension); //initial guess types::RMatrix natural_guess(1, dimension, 0.); // realspace coordinates provided by initial guess - types::RMatrix physical_guess(1,dimension); + types::RMatrix physical_guess(1, dimension); // objective function f = real_coords - physical_guess types::RMatrix f(1,dimension); // dnds computed on the natural_guess types::RMatrix dnds(nb_nodes_per_element,spatial_dimension); // dxds computed on the natural_guess types::RMatrix dxds(spatial_dimension,dimension); // transposed dxds computed on the natural_guess types::RMatrix dxds_t(dimension,spatial_dimension); // G = dxds * dxds_t types::RMatrix G(spatial_dimension,spatial_dimension); // Ginv = G^{-1} types::RMatrix Ginv(spatial_dimension,spatial_dimension); // J = Ginv * dxds types::RMatrix J(spatial_dimension,dimension); // dxi = \xi_{k+1} - \xi in the iterative process types::RMatrix dxi(1,spatial_dimension); /* --------------------------- */ // init before iteration loop /* --------------------------- */ // do interpolation interpolateOnNaturalCoordinates(natural_guess.storage(),node_coords.storage(),dimension,physical_guess.storage()); // compute initial objective function value f = real_coords - physical_guess f = physical_guess; f*= -1.; f+= mreal_coords; // compute initial error Real inverse_map_error = f.norm(); /* --------------------------- */ // iteration loop /* --------------------------- */ while(tolerance < inverse_map_error) { //compute dxds computeDNDS(natural_guess.storage(), dnds.storage()); computeDXDS(dnds.storage(),node_coords.storage(),dimension,dxds.storage()); //compute G dxds_t = dxds; dxds_t.transpose(); G.mul(dxds,dxds_t); // inverse G if (spatial_dimension == 1) Ginv[0] = 1./G[0]; else if (spatial_dimension == 2) Math::inv2(G.storage(),Ginv.storage()); else if (spatial_dimension == 3) Math::inv3(G.storage(),Ginv.storage()); //compute J J.mul(Ginv,dxds); //compute increment dxi.mul(f,J); //update our guess natural_guess += dxi; //interpolate interpolateOnNaturalCoordinates(natural_guess.storage(),node_coords.storage(),dimension,physical_guess.storage()); // compute error f = physical_guess; f*= -1.; f+= mreal_coords; inverse_map_error = f.norm(); } memcpy(natural_coords.storage(), natural_guess.storage(), sizeof(Real) * natural_coords.size()); } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapes(const Real * natural_coords, Real * shapes) { computeShapes(natural_coords, shapes, NULL,0); } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapes(__attribute__ ((unused)) const Real * natural_coords, __attribute__ ((unused)) Real * shapes, __attribute__ ((unused)) const Real * local_coord, __attribute__ ((unused)) UInt id) { AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeDNDS(__attribute__((unused)) const Real * natural_coords, __attribute__((unused)) Real * dnds){ // computeDNDS(natural_coords, dnds); AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapeDerivatives(__attribute__ ((unused)) const Real * natural_coords, __attribute__ ((unused)) Real * shape_deriv, __attribute__ ((unused)) const Real * local_coord, __attribute__ ((unused)) UInt id) { AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeJacobian(__attribute__ ((unused)) const Real * dxds, __attribute__ ((unused)) const UInt dimension, __attribute__ ((unused)) Real & jac) { //AKANTU_DEBUG_ERROR("Function not implemented for type : " << type); } /* -------------------------------------------------------------------------- */ template inline Real * ElementClass::getGaussIntegrationWeights() { return gauss_integration_weights; } /* -------------------------------------------------------------------------- */ template inline bool ElementClass::contains(__attribute__ ((unused)) const types::RVector & natural_coords) { AKANTU_DEBUG_TO_IMPLEMENT(); return false; } /* -------------------------------------------------------------------------- */ #include "element_classes/element_class_segment_2_inline_impl.cc" #include "element_classes/element_class_segment_3_inline_impl.cc" #include "element_classes/element_class_triangle_3_inline_impl.cc" #include "element_classes/element_class_triangle_6_inline_impl.cc" #include "element_classes/element_class_tetrahedron_4_inline_impl.cc" #include "element_classes/element_class_tetrahedron_10_inline_impl.cc" #include "element_classes/element_class_quadrangle_4_inline_impl.cc" #include "element_classes/element_class_quadrangle_8_inline_impl.cc" #include "element_classes/element_class_hexahedron_8_inline_impl.cc" #include "element_classes/element_class_bernoulli_beam_2_inline_impl.cc" diff --git a/src/fem/mesh.hh b/src/fem/mesh.hh index b47b6b8bd..b37520257 100644 --- a/src/fem/mesh.hh +++ b/src/fem/mesh.hh @@ -1,595 +1,596 @@ /** * @file mesh.hh * @author Nicolas Richart * @date Wed Jun 16 11:53:53 2010 * * @brief the class representing the meshes * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_HH__ #define __AKANTU_MESH_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "aka_vector.hh" #include "element_class.hh" #include "by_element_type.hh" #include "aka_event_handler.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* Element */ /* -------------------------------------------------------------------------- */ class Element { public: Element(ElementType type = _not_defined, UInt element = 0, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) : type(type), element(element), ghost_type(ghost_type), kind(kind) {}; Element(const Element & element) { this->type = element.type; this->element = element.element; this->ghost_type = element.ghost_type; this->kind = element.kind; } inline bool operator==(const Element & elem) const { return ((element == elem.element) && (type == elem.type) && (ghost_type == elem.ghost_type) && (kind == elem.kind)); } inline bool operator!=(const Element & elem) const { return ((element != elem.element) || (type != elem.type) || (ghost_type != elem.ghost_type) || (kind != elem.kind)); } bool operator<(const Element& rhs) const { bool res = ((this->kind < rhs.kind) || ((this->kind == rhs.kind) && ((this->ghost_type < rhs.ghost_type) || ((this->ghost_type == rhs.ghost_type) && ((this->type < rhs.type) || ((this->type == rhs.type) && (this->element < rhs.element))))))); return res; } virtual ~Element() {}; /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; public: ElementType type; UInt element; GhostType ghost_type; ElementKind kind; }; struct CompElementLess { bool operator() (const Element& lhs, const Element& rhs) const { return lhs < rhs; } }; extern const Element ElementNull; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* Mesh modifications events */ /* -------------------------------------------------------------------------- */ template class MeshEvent { public: const Vector & getList() const { return list; } Vector & getList() { return list; } protected: Vector list; }; class Mesh; class NewNodesEvent : public MeshEvent { }; class RemovedNodesEvent : public MeshEvent { public: inline RemovedNodesEvent(const Mesh & mesh); AKANTU_GET_MACRO_NOT_CONST(NewNumbering, new_numbering, Vector &); AKANTU_GET_MACRO(NewNumbering, new_numbering, const Vector &); private: Vector new_numbering; }; class NewElementsEvent : public MeshEvent { }; class RemovedElementsEvent : public MeshEvent { public: inline RemovedElementsEvent(const Mesh & mesh); AKANTU_GET_MACRO(NewNumbering, new_numbering, const ByElementTypeUInt &); AKANTU_GET_MACRO_NOT_CONST(NewNumbering, new_numbering, ByElementTypeUInt &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(NewNumbering, new_numbering, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(NewNumbering, new_numbering, UInt); protected: ByElementTypeUInt new_numbering; }; /* -------------------------------------------------------------------------- */ class MeshEventHandler { public: virtual ~MeshEventHandler() {}; /* ------------------------------------------------------------------------ */ /* Internal code */ /* ------------------------------------------------------------------------ */ public: inline void sendEvent(const NewNodesEvent & event) { onNodesAdded (event.getList()); } inline void sendEvent(const RemovedNodesEvent & event) { onNodesRemoved(event.getList(), event.getNewNumbering()); } inline void sendEvent(const NewElementsEvent & event) { onElementsAdded (event.getList()); } inline void sendEvent(const RemovedElementsEvent & event) { onElementsRemoved(event.getList(), event.getNewNumbering()); } /* ------------------------------------------------------------------------ */ /* Interface */ /* ------------------------------------------------------------------------ */ public: virtual void onNodesAdded (__attribute__((unused)) const Vector & nodes_list) { } virtual void onNodesRemoved(__attribute__((unused)) const Vector & nodes_list, __attribute__((unused)) const Vector & new_numbering) { } virtual void onElementsAdded (__attribute__((unused)) const Vector & elements_list) { } virtual void onElementsRemoved(__attribute__((unused)) const Vector & elements_list, __attribute__((unused)) const ByElementTypeUInt & new_numbering) { } }; /* -------------------------------------------------------------------------- */ /* Mesh */ /* -------------------------------------------------------------------------- */ /** * @class Mesh this contain the coordinates of the nodes in the Mesh.nodes * Vector, and the connectivity. The connectivity are stored in by element * types. * * To know all the element types present in a mesh you can get the * Mesh::ConnectivityTypeList * * In order to loop on all element you have to loop on all types like this : * @code Mesh::type_iterator it = mesh.firstType(dim, ghost_type); Mesh::type_iterator end = mesh.lastType(dim, ghost_type); for(; it != end; ++it) { UInt nb_element = mesh.getNbElement(*it); const Vector & conn = mesh.getConnectivity(*it); for(UInt e = 0; e < nb_element; ++e) { ... } } @endcode */ class Mesh : protected Memory, public EventHandlerManager { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /// constructor that create nodes coordinates array Mesh(UInt spatial_dimension, const ID & id = "mesh", const MemoryID & memory_id = 0); /// constructor that use an existing nodes coordinates array, by knowing its ID Mesh(UInt spatial_dimension, const ID & nodes_id, const ID & id = "mesh", const MemoryID & memory_id = 0); /** * constructor that use an existing nodes coordinates * array, by getting the vector of coordinates */ Mesh(UInt spatial_dimension, Vector & nodes, const ID & id = "mesh", const MemoryID & memory_id = 0); virtual ~Mesh(); /// @typedef ConnectivityTypeList list of the types present in a Mesh typedef std::set ConnectivityTypeList; // typedef Vector * NormalsMap[_max_element_type]; private: /// initialize the connectivity to NULL and other stuff void init(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; /// function that computes the bounding box (fills xmin, xmax) void computeBoundingBox(); /// init a by-element-type real vector with provided ids template void initByElementTypeVector(ByElementTypeVector & v, UInt nb_component, UInt spatial_dimension, const bool & flag_nb_node_per_elem_multiply = false, ElementKind element_kind = _ek_regular, bool size_to_nb_element = false) const; /// @todo: think about nicer way to do it /// extract coordinates of nodes from an element template inline void extractNodalValuesFromElement(const Vector & nodal_values, T * elemental_values, UInt * connectivity, UInt n_nodes, UInt nb_degree_of_freedom) const; /// extract coordinates of nodes from a reversed element inline void extractNodalCoordinatesFromPBCElement(Real * local_coords, UInt * connectivity, UInt n_nodes); /// convert a element to a linearized element inline UInt elementToLinearized(const Element & elem); /// convert a linearized element to an element inline Element linearizedToElement (UInt linearized_element); /// update the types offsets array for the conversions inline void updateTypesOffsets(const GhostType & ghost_type); /// add a Vector of connectivity for the type . inline void addConnectivityType(const ElementType & type); /* ------------------------------------------------------------------------ */ template inline void sendEvent(Event & event) { if(event.getList().getSize() != 0) EventHandlerManager::sendEvent(event); } /* ------------------------------------------------------------------------ */ template inline void removeNodesFromVector(Vector & vect, const Vector & new_numbering); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(ID, id, const ID &); /// get the spatial dimension of the mesh = number of component of the coordinates AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the nodes Vector aka coordinates AKANTU_GET_MACRO(Nodes, *nodes, const Vector &); AKANTU_GET_MACRO_NOT_CONST(Nodes, *nodes, Vector &); /// get the number of nodes AKANTU_GET_MACRO(NbNodes, nodes->getSize(), UInt); /// get the Vector of global ids of the nodes (only used in parallel) AKANTU_GET_MACRO(GlobalNodesIds, *nodes_global_ids, const Vector &); /// get the global id of a node inline UInt getNodeGlobalId(UInt local_id) const; /// get the global number of nodes inline UInt getNbGlobalNodes() const; /// get the nodes type Vector AKANTU_GET_MACRO(NodesType, *nodes_type, const Vector &); inline Int getNodeType(UInt local_id) const; /// say if a node is a pure ghost node inline bool isPureGhostNode(UInt n) const; /// say if a node is pur local or master node inline bool isLocalOrMasterNode(UInt n) const; inline bool isLocalNode(UInt n) const; inline bool isMasterNode(UInt n) const; inline bool isSlaveNode(UInt n) const; AKANTU_GET_MACRO(XMin, lower_bounds[0], Real); AKANTU_GET_MACRO(YMin, lower_bounds[1], Real); AKANTU_GET_MACRO(ZMin, lower_bounds[2], Real); AKANTU_GET_MACRO(XMax, upper_bounds[0], Real); AKANTU_GET_MACRO(YMax, upper_bounds[1], Real); AKANTU_GET_MACRO(ZMax, upper_bounds[2], Real); inline void getLowerBounds(Real * lower) const; inline void getUpperBounds(Real * upper) const; inline void getLocalLowerBounds(Real * lower) const; inline void getLocalUpperBounds(Real * upper) const; /// get the number of surfaces AKANTU_GET_MACRO(NbSurfaces, nb_surfaces, UInt); /// get the connectivity Vector for a given type AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Connectivity, connectivities, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Connectivity, connectivities, UInt); + AKANTU_GET_MACRO(Connectivities, connectivities, const ByElementTypeVector &); /// @todo take out this set, if mesh can read surface id /// set the number of surfaces AKANTU_SET_MACRO(NbSurfaces, nb_surfaces, UInt); /// get the number of element of a type in the mesh inline UInt getNbElement(const ElementType & type, const GhostType & ghost_type = _not_ghost) const; // /// get the number of ghost element of a type in the mesh // inline UInt getNbGhostElement(const ElementType & type) const; /// get the connectivity list either for the elements or the ghost elements inline const ConnectivityTypeList & getConnectivityTypeList(const GhostType & ghost_type = _not_ghost) const; /// compute the barycenter of a given element inline void getBarycenter(UInt element, const ElementType & type, Real * barycenter, GhostType ghost_type = _not_ghost) const; /// get the element connected to a subelement AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementToSubelement, element_to_subelement, Vector); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementToSubelement, element_to_subelement, Vector); /// get the subelement connected to an element AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(SubelementToElement, subelement_to_element, Element); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(SubelementToElement, subelement_to_element, Element); /// get the surface values of facets AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(SurfaceID, surface_id, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(SurfaceID, surface_id, UInt); /// set the int data to the surface id vectors void setSurfaceIDsFromIntData(const std::string & data_name); inline const Vector & getUIntData(const ElementType & el_type, const std::string & data_name, const GhostType & ghost_type = _not_ghost) const; /* ------------------------------------------------------------------------ */ /* Wrappers on ElementClass functions */ /* ------------------------------------------------------------------------ */ public: /// get the number of nodes per element for a given element type static inline UInt getNbNodesPerElement(const ElementType & type); /// get the number of nodes per element for a given element type considered as /// a first order element static inline ElementType getP1ElementType(const ElementType & type); /// get the kind of the element type static inline ElementKind getKind(const ElementType & type); /// get spatial dimension of a type of element static inline UInt getSpatialDimension(const ElementType & type); /// get number of facets of a given element type static inline UInt getNbFacetsPerElement(const ElementType & type); /// get local connectivity of a facet for a given facet type static inline UInt ** getFacetLocalConnectivity(const ElementType & type); /// get the type of the surface element associated to a given element static inline ElementType getFacetElementType(const ElementType & type); /// get the pointer to the list of elements for a given type inline Vector * getReversedElementsPBCPointer(const ElementType & type); /* ------------------------------------------------------------------------ */ /* Element type Iterator */ /* ------------------------------------------------------------------------ */ typedef ByElementTypeUInt::type_iterator type_iterator; inline type_iterator firstType(UInt dim = 0, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const { return connectivities.firstType(dim, ghost_type, kind); } inline type_iterator lastType(UInt dim = 0, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const { return connectivities.lastType(dim, ghost_type, kind); } /* ------------------------------------------------------------------------ */ /* Private methods for friends */ /* ------------------------------------------------------------------------ */ private: friend class MeshIOMSH; friend class MeshIOMSHStruct; friend class MeshIODiana; friend class MeshUtils; friend class DistributedSynchronizer; AKANTU_GET_MACRO(NodesPointer, nodes, Vector *); /// get a pointer to the nodes_global_ids Vector and create it if necessary inline Vector * getNodesGlobalIdsPointer(); /// get a pointer to the nodes_type Vector and create it if necessary inline Vector * getNodesTypePointer(); /// get a pointer to the connectivity Vector for the given type and create it if necessary inline Vector * getConnectivityPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); // inline Vector * getNormalsPointer(ElementType type) const; /// get a pointer to the surface_id Vector for the given type and create it if necessary inline Vector * getSurfaceIDPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); /// get the UIntDataMap for a given ElementType inline UIntDataMap & getUIntDataMap(const ElementType & el_type, const GhostType & ghost_type = _not_ghost); /// get the IntDataMap pointer (modifyable) for a given ElementType inline Vector * getUIntDataPointer(const ElementType & el_type, const std::string & data_name, const GhostType & ghost_type = _not_ghost); /// get a pointer to the element_to_subelement Vector for the given type and create it if necessary inline Vector > * getElementToSubelementPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); /// get a pointer to the subelement_to_element Vector for the given type and create it if necessary inline Vector * getSubelementToElementPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// id of the mesh ID id; /// array of the nodes coordinates Vector * nodes; /// global node ids Vector * nodes_global_ids; /// node type, -3 pure ghost, -2 master for the node, -1 normal node, i in /// [0-N] slave node and master is proc i Vector * nodes_type; /// global number of nodes; UInt nb_global_nodes; /// boolean to know if the nodes have to be deleted with the mesh or not bool created_nodes; /// all class of elements present in this mesh (for heterogenous meshes) ByElementTypeUInt connectivities; /// map to normals for all class of elements present in this mesh ByElementTypeReal normals; /// list of all existing types in the mesh ConnectivityTypeList type_set; /// the spatial dimension of this mesh UInt spatial_dimension; /// types offsets Vector types_offsets; /// list of all existing types in the mesh ConnectivityTypeList ghost_type_set; /// ghost types offsets Vector ghost_types_offsets; /// number of surfaces present in this mesh UInt nb_surfaces; /// surface id of the surface elements in this mesh ByElementTypeUInt surface_id; /// min of coordinates Real lower_bounds[3]; /// max of coordinates Real upper_bounds[3]; /// size covered by the mesh on each direction Real size[3]; /// local min of coordinates Real local_lower_bounds[3]; /// local max of coordinates Real local_upper_bounds[3]; /// List of elements connected to subelements ByElementTypeVector > element_to_subelement; /// List of subelements connected to elements ByElementTypeVector subelement_to_element; // /// list of elements that are reversed due to pbc // ByElementTypeUInt reversed_elements_pbc; // /// direction in which pbc are to be applied // UInt pbc_directions[3]; /// list of the vectors corresponding to tags in the mesh ByElementTypeUIntDataMap uint_data; }; /* -------------------------------------------------------------------------- */ /* Inline functions */ /* -------------------------------------------------------------------------- */ /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Element & _this) { _this.printself(stream); return stream; } #if defined (AKANTU_INCLUDE_INLINE_IMPL) # include "mesh_inline_impl.cc" #endif #include "by_element_type_tmpl.hh" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Mesh & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_MESH_HH__ */ diff --git a/src/mesh_utils/mesh_utils.cc b/src/mesh_utils/mesh_utils.cc index de1f5884f..4ed00d518 100644 --- a/src/mesh_utils/mesh_utils.cc +++ b/src/mesh_utils/mesh_utils.cc @@ -1,919 +1,915 @@ /** * @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 * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "mesh_utils.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ void MeshUtils::buildNode2Elements(const Mesh & mesh, CSR & node_to_elem, UInt spatial_dimension) { AKANTU_DEBUG_IN(); if (spatial_dimension == 0) spatial_dimension = mesh.getSpatialDimension(); UInt nb_nodes = mesh.getNbNodes(); const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; UInt nb_types = type_list.size(); UInt nb_good_types = 0; UInt nb_nodes_per_element[nb_types]; // UInt nb_nodes_per_element_p1[nb_types]; UInt * conn_val[nb_types]; UInt nb_element[nb_types]; // ElementType type_p1; for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; if(Mesh::getSpatialDimension(type) != spatial_dimension) continue; nb_nodes_per_element[nb_good_types] = Mesh::getNbNodesPerElement(type); // type_p1 = Mesh::getP1ElementType(type); // nb_nodes_per_element_p1[nb_good_types] = Mesh::getNbNodesPerElement(type_p1); conn_val[nb_good_types] = mesh.getConnectivity(type, _not_ghost).values; nb_element[nb_good_types] = mesh.getConnectivity(type, _not_ghost).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_to_elem.resizeRows(nb_nodes); node_to_elem.clearRows(); // node_offset.resize(nb_nodes + 1); // UInt * node_offset_val = node_offset.values; /// count number of occurrence of each node // memset(node_offset_val, 0, (nb_nodes + 1)*sizeof(UInt)); for (UInt t = 0; t < nb_good_types; ++t) { for (UInt el = 0; el < nb_element[t]; ++el) { UInt el_offset = el*nb_nodes_per_element[t]; for (UInt n = 0; n < nb_nodes_per_element[t]; ++n) { ++node_to_elem.rowOffset(conn_val[t][el_offset + 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; node_to_elem.countToCSR(); node_to_elem.resizeCols(); node_to_elem.beginInsertions(); /// 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[t]; ++n) node_to_elem.insertInRow(conn_val[t][el_offset + n], linearized_el); // node_to_elem_val[node_offset_val[conn_val[t][el_offset + n]]++] = linearized_el; } node_to_elem.endInsertions(); // for (UInt i = nb_nodes; i > 0; --i) node_offset_val[i] = node_offset_val[i-1]; // node_offset_val[0] = 0; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildNode2ElementsByElementType(const Mesh & mesh, ElementType type, CSR & node_to_elem) { // Vector & node_offset, // Vector & node_to_elem) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_elements = mesh.getConnectivity(type, _not_ghost).getSize(); UInt * conn_val = mesh.getConnectivity(type, _not_ghost).values; /// array for the node-element list node_to_elem.resizeRows(nb_nodes); node_to_elem.clearRows(); // node_offset.resize(nb_nodes + 1); // UInt * node_offset_val = node_offset.values; // memset(node_offset_val, 0, (nb_nodes + 1)*sizeof(UInt)); /// count number of occurrence of each node for (UInt el = 0; el < nb_elements; ++el) { UInt el_offset = el*nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) ++node_to_elem.rowOffset(conn_val[el_offset + n]); // node_offset_val[conn_val[nb_nodes_per_element*el + n]]++; } /// convert the occurrence array in a csr one // for (UInt i = 1; i < nb_nodes; ++i) node_offset_val[i] += node_offset_val[i-1]; // for (UInt i = nb_nodes; i > 0; --i) node_offset_val[i] = node_offset_val[i-1]; // node_offset_val[0] = 0; node_to_elem.countToCSR(); node_to_elem.resizeCols(); node_to_elem.beginInsertions(); /// save the element index in the node-element list // node_to_elem.resize(node_offset_val[nb_nodes]); // UInt * node_to_elem_val = node_to_elem.values; for (UInt el = 0; el < nb_elements; ++el) { UInt el_offset = el*nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) { node_to_elem.insertInRow(conn_val[el_offset + n], el); } // node_to_elem_val[node_offset_val[conn_val[nb_nodes_per_element*el + n]]] = el; // node_offset_val[conn_val[nb_nodes_per_element*el + n]]++; } // /// rearrange node_offset to start with 0 // for (UInt i = nb_nodes; i > 0; --i) node_offset_val[i] = node_offset_val[i-1]; // node_offset_val[0] = 0; node_to_elem.endInsertions(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildFacets(Mesh & mesh){ AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); ByElementTypeReal barycenter; buildFacetsDimension(mesh, mesh, true, spatial_dimension, barycenter); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildAllFacets(Mesh & mesh, Mesh & mesh_facets) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); ByElementTypeReal barycenter; /// generate facets buildFacetsDimension(mesh, mesh_facets, false, spatial_dimension, barycenter); /// compute their barycenters mesh_facets.initByElementTypeVector(barycenter, spatial_dimension, spatial_dimension - 1); Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1); Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1); for(; it != end; ++it) { UInt nb_element = mesh_facets.getNbElement(*it); barycenter(*it).resize(nb_element); Vector::iterator bary = barycenter(*it).begin(spatial_dimension); Vector::iterator bary_end = barycenter(*it).end(spatial_dimension); for (UInt el = 0; bary != bary_end; ++bary, ++el) { mesh_facets.getBarycenter(el, *it, bary->storage()); } } /// sort facets and generate subfacets for (UInt i = spatial_dimension - 1; i > 0; --i) { buildFacetsDimension(mesh_facets, mesh_facets, false, i, barycenter); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildFacetsDimension(Mesh & mesh, Mesh & mesh_facets, bool boundary_only, UInt dimension, ByElementTypeReal & barycenter){ AKANTU_DEBUG_IN(); const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; UInt nb_types = type_list.size(); UInt nb_good_types = 0; UInt spatial_dimension = mesh.getSpatialDimension(); 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 > * element_to_subelement[nb_types]; Vector * subelement_to_element[nb_types]; UInt * conn_val[nb_types]; UInt nb_element[nb_types]; ElementType facet_type; Real epsilon = std::numeric_limits::epsilon(); for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; if(mesh.getSpatialDimension(type) != dimension) continue; nb_nodes_per_element[nb_good_types] = mesh.getNbNodesPerElement(type); facet_type = mesh.getFacetElementType(type); nb_facets[nb_good_types] = mesh.getNbFacetsPerElement(type); node_in_facet[nb_good_types] = mesh.getFacetLocalConnectivity(type); nb_nodes_per_facet[nb_good_types] = mesh.getNbNodesPerElement(facet_type); // getting connectivity of boundary facets connectivity_facets[nb_good_types] = mesh_facets.getConnectivityPointer(facet_type); connectivity_facets[nb_good_types]->resize(0); element_to_subelement[nb_good_types] = mesh_facets.getElementToSubelementPointer(facet_type); element_to_subelement[nb_good_types]->resize(0); conn_val[nb_good_types] = mesh.getConnectivity(type, _not_ghost).values; nb_element[nb_good_types] = mesh.getConnectivity(type, _not_ghost).getSize(); subelement_to_element[nb_good_types] = mesh_facets.getSubelementToElementPointer(type); subelement_to_element[nb_good_types]->resize(nb_element[nb_good_types]); nb_good_types++; } CSR node_to_elem; // Vector node_offset; // Vector node_to_elem; // buildNode2Elements(mesh,node_offset,node_to_elem); buildNode2Elements(mesh, node_to_elem, dimension); // std::cout << "node offset " << std::endl << node_offset << std::endl; // std::cout << "node to elem " << std::endl << node_to_elem << std::endl; 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 CSR::iterator first_node_elements; //UInt * first_node_elements = node_to_elem.values+node_offset.values[facet_nodes[0]]; UInt first_node_nelements = node_to_elem.getNbCols(facet_nodes[0]); // node_offset.values[facet_nodes[0]+1]- // node_offset.values[facet_nodes[0]]; counter.resize(first_node_nelements); counter.clear(); //loop over the other nodes to search intersecting elements, //which are the elements that share another node with the //starting element after first_node CSR::iterator first_node_elements_end; first_node_elements_end = node_to_elem.end(facet_nodes[0]); for (UInt n = 1; n < nb_nodes_per_facet[t]; ++n) { CSR::iterator node_elements, node_elements_begin, node_elements_end; node_elements_begin = node_to_elem.begin(facet_nodes[n]); node_elements_end = node_to_elem.end(facet_nodes[n]); // UInt * node_elements = node_to_elem.values+node_offset.values[facet_nodes[n]]; // node_offset.values[facet_nodes[n]+1]- // node_offset.values[facet_nodes[n]]; UInt local_el = 0; for (first_node_elements = node_to_elem.begin(facet_nodes[0]); first_node_elements != first_node_elements_end; ++first_node_elements, ++local_el) { for (node_elements = node_elements_begin; node_elements != node_elements_end; ++node_elements) { if (*first_node_elements == *node_elements) { ++counter.values[local_el]; // it may cause trouble: break; } } // if (counter.values[local_el] == 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; // counting the number of elements connected to the facets and // taking the minimum element number, because the facet should // be inserted just once UInt nb_element_connected_to_facet = 0; UInt minimum_el_index = std::numeric_limits::max(); Vector connected_elements; for (UInt el1 = 0; el1 < counter.getSize(); el1++) { UInt el_index = node_to_elem(facet_nodes[0], el1); if (counter.values[el1] == nb_nodes_per_facet[t]-1) { ++nb_element_connected_to_facet; minimum_el_index = std::min(minimum_el_index, el_index); connected_elements.push_back(el_index); } } if (minimum_el_index == linearized_el) { if (!boundary_only || (boundary_only && nb_element_connected_to_facet == 1)) { connectivity_facets[t]->push_back(facet_nodes); UInt current_nb_facets = element_to_subelement[t]->getSize(); element_to_subelement[t]->resize(current_nb_facets + 1); Vector & elements = (*element_to_subelement[t])(current_nb_facets); // build elements_on_facets: linearized_el must come first // in order to store the facet in the correct direction // and avoid to invert the sign in the normal computation elements.push_back(mesh.linearizedToElement(linearized_el)); /// boundary facet if (nb_element_connected_to_facet == 1) elements.push_back(ElementNull); /// internal facet else if (nb_element_connected_to_facet == 2) elements.push_back(mesh.linearizedToElement(connected_elements.values[1])); /// facet of facet else { UInt nb_connected = connected_elements.getSize(); for (UInt i = 1; i < nb_connected; ++i) { elements.push_back(mesh.linearizedToElement(connected_elements.values[i])); } /// check if sorting is needed: /// - in 3D to sort triangles around segments /// - in 2D to sort segments around points if (dimension == spatial_dimension - 1) { /// barycentrical coordinates for each connected /// element with respect to start_node Vector connected_nodes(nb_connected, 2); const Vector & coord = mesh_facets.getNodes(); const Vector & facet_conn = mesh_facets.getConnectivity(facet_type); /// node around which the sorting is carried out is /// the first node of the current facet UInt start_node = facet_conn(facet_conn.getSize()-1, 0); Real start_coord[spatial_dimension]; for (UInt dim = 0; dim < spatial_dimension; ++dim) { start_coord[dim] = coord(start_node, dim); } if (spatial_dimension == 3) { /// vector connecting facet first node to second Real tangent[spatial_dimension]; /// vector connecting facet first node and /// barycenter of elements(0) Real temp[spatial_dimension]; /// two normals to the segment facet to define the /// reference system Real normal1[spatial_dimension]; Real normal2[spatial_dimension]; Vector & bar = barycenter(elements(0).type); /// facet second node UInt second_node = facet_conn(facet_conn.getSize()-1, 1); /// construction of tangent and temp arrays for (UInt dim = 0; dim < spatial_dimension; ++dim) { Real x1, x2; x1 = coord(second_node, dim); x2 = bar(elements(0).element, dim); tangent[dim] = x1 - start_coord[dim]; temp[dim] = x2 - start_coord[dim]; } /// get normal1 and normal2 Math::normalize3(tangent); Math::vectorProduct3(tangent, temp, normal1); Math::normalize3(normal1); Math::vectorProduct3(tangent, normal1, normal2); for (UInt n = 0; n < nb_connected; ++n) { Real bary_coord[spatial_dimension]; Vector & bary = barycenter(elements(n).type); /// get the barycenter local coordinates for (UInt dim = 0; dim < spatial_dimension; ++dim) { bary_coord[dim] = bary(elements(n).element, dim) - start_coord[dim]; } /// project the barycenter coordinates on the two /// normals to have them on the same plane connected_nodes(n, 0) = Math::vectorDot(bary_coord, normal1, spatial_dimension); connected_nodes(n, 1) = Math::vectorDot(bary_coord, normal2, spatial_dimension); } } else if (spatial_dimension == 2) { for (UInt n = 0; n < nb_connected; ++n) { Vector & bary = barycenter(elements(n).type); /// get the barycenter local coordinates for (UInt dim = 0; dim < spatial_dimension; ++dim) { connected_nodes(n, dim) = bary(elements(n).element, dim) - start_coord[dim]; } } } /// associate to each element a real value based on /// atan2 function (check wikipedia) std::map atan2; for (UInt n = 0; n < nb_connected; ++n) { Real x = connected_nodes(n, 0); Real y = connected_nodes(n, 1); /// in order to avoid division by zero: if (std::abs(y) <= std::abs(y) * epsilon && x < 0) y = Math::getTolerance(); atan2[elements(n)] = y / (sqrt(x * x + y * y) + x); } /// sort elements according to their atan2 values ElementSorter sorter(atan2); std::sort(elements.storage(), elements.storage() + elements.getSize(), sorter); } } /// current facet index UInt current_facet = connectivity_facets[t]->getSize() - 1; /// loop on every element connected to current facet and /// insert current facet in the first free spot of the /// subelement_to_element vector for (UInt elem = 0; elem < elements.getSize(); ++elem) { if (elements(elem).type != _not_defined) { for (UInt f_in = 0; f_in < nb_facets[t]; ++f_in) { if ((*subelement_to_element[t])(elements(elem).element, f_in).type == _not_defined) { (*subelement_to_element[t])(elements(elem).element, f_in).type = facet_type; (*subelement_to_element[t])(elements(elem).element, f_in).element = current_facet; break; } } } } } } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ // void MeshUtils::buildNormals(Mesh & mesh,UInt spatial_dimension){ // AKANTU_DEBUG_IN(); // const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); // Mesh::ConnectivityTypeList::const_iterator it; // UInt nb_types = type_list.size(); // UInt nb_nodes_per_element[nb_types]; // UInt nb_nodes_per_element_p1[nb_types]; // UInt nb_good_types = 0; // Vector * connectivity[nb_types]; // Vector * normals[nb_types]; // if (spatial_dimension == 0) spatial_dimension = mesh.getSpatialDimension(); // for(it = type_list.begin(); it != type_list.end(); ++it) { // ElementType type = *it; // ElementType type_p1 = mesh.getP1ElementType(type); // if(mesh.getSpatialDimension(type) != spatial_dimension) continue; // nb_nodes_per_element[nb_good_types] = mesh.getNbNodesPerElement(type); // nb_nodes_per_element_p1[nb_good_types] = mesh.getNbNodesPerElement(type_p1); // // getting connectivity // connectivity[nb_good_types] = mesh.getConnectivityPointer(type); // if (!connectivity[nb_good_types]) // AKANTU_DEBUG_ERROR("connectivity is not allocatted : this should probably not have happened"); // //getting array of normals // normals[nb_good_types] = mesh.getNormalsPointer(type); // if(normals[nb_good_types]) // normals[nb_good_types]->resize(0); // else // normals[nb_good_types] = &mesh.createNormals(type); // nb_good_types++; // } // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ void MeshUtils::renumberMeshNodes(Mesh & mesh, UInt * local_connectivities, UInt nb_local_element, UInt nb_ghost_element, ElementType type, - Vector & nodes_numbers) { + Vector & old_nodes_numbers) { AKANTU_DEBUG_IN(); - nodes_numbers.resize(0); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); std::map renumbering_map; + old_nodes_numbers.resize(0); + /// renumber the nodes renumberNodesInConnectivity(local_connectivities, (nb_local_element + nb_ghost_element)*nb_nodes_per_element, renumbering_map, - nodes_numbers); + old_nodes_numbers); renumbering_map.clear(); /// copy the renumbered connectivity to the right place Vector * local_conn = mesh.getConnectivityPointer(type); local_conn->resize(nb_local_element); memcpy(local_conn->values, local_connectivities, nb_local_element * nb_nodes_per_element * sizeof(UInt)); Vector * ghost_conn = mesh.getConnectivityPointer(type,_ghost); ghost_conn->resize(nb_ghost_element); memcpy(ghost_conn->values, local_connectivities + nb_local_element * nb_nodes_per_element, nb_ghost_element * nb_nodes_per_element * sizeof(UInt)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::renumberNodesInConnectivity(UInt * list_nodes, UInt nb_nodes, std::map & renumbering_map, Vector & nodes_numbers) { AKANTU_DEBUG_IN(); UInt * connectivity = list_nodes; UInt new_node_num = renumbering_map.size(); for (UInt n = 0; n < nb_nodes; ++n, ++connectivity) { UInt & node = *connectivity; std::map::iterator it = renumbering_map.find(node); if(it == renumbering_map.end()) { UInt old_node = node; - nodes_numbers(node) = new_node_num; - node = new_node_num; + nodes_numbers.push_back(old_node); renumbering_map[old_node] = new_node_num; + node = new_node_num; ++new_node_num; } else { node = it->second; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::purifyMesh(Mesh & mesh) { AKANTU_DEBUG_IN(); std::map renumbering_map; RemovedNodesEvent remove_nodes(mesh); Vector & nodes_removed = remove_nodes.getList(); - Vector & node_numbers = remove_nodes.getNewNumbering(); - std::fill(node_numbers.begin(), node_numbers.end(), UInt(-1)); + Vector node_numbers; for (UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type = (GhostType) gt; Mesh::type_iterator it = mesh.firstType(0, ghost_type, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(0, ghost_type, _ek_not_defined); for(; it != end; ++it) { ElementType type(*it); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const Vector & connectivity_vect = mesh.getConnectivity(type, ghost_type); UInt nb_element(connectivity_vect.getSize()); UInt * connectivity = connectivity_vect.storage(); renumberNodesInConnectivity (connectivity, nb_element*nb_nodes_per_element, renumbering_map, node_numbers); } } - for (UInt i = 0; i < node_numbers.getSize(); ++i) { - if(node_numbers(i) == UInt(-1)) - nodes_removed.push_back(i); - } + Vector & new_numbering = remove_nodes.getNewNumbering(); + std::fill(new_numbering.begin(), new_numbering.end(), UInt(-1)); - mesh.sendEvent(remove_nodes); - - // UInt spatial_dimension = mesh.getSpatialDimension(); - // Vector & nodes = const_cast &>(mesh.getNodes()); - // Vector::iterator nodes_it = nodes.begin(spatial_dimension); - // Vector new_nodes(0, spatial_dimension); + std::map::iterator it = renumbering_map.begin(); + std::map::iterator end = renumbering_map.end(); + for (; it != end; ++it) { + new_numbering(it->first) = it->second; + } - // for (UInt i = 0; i < node_numbers.getSize(); ++i) { - // new_nodes.push_back(nodes_it + node_numbers(i)); - // } - // std::copy(new_nodes.begin(spatial_dimension), new_nodes.end(spatial_dimension), - // nodes.begin(spatial_dimension)); + for (UInt i = 0; i < new_numbering.getSize(); ++i) { + if(new_numbering(i) == UInt(-1)) + nodes_removed.push_back(i); + } - // nodes.resize(node_numbers.getSize()); + mesh.sendEvent(remove_nodes); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::setUIntData(Mesh & mesh, UInt * data, UInt nb_tags, const ElementType & type) { AKANTU_DEBUG_IN(); UInt nb_element = mesh.getNbElement(type, _not_ghost); UInt nb_ghost_element = mesh.getNbElement(type, _ghost); char * names = reinterpret_cast(data + (nb_element + nb_ghost_element) * nb_tags); UIntDataMap & uint_data_map = mesh.getUIntDataMap(type, _not_ghost); UIntDataMap & ghost_uint_data_map = mesh.getUIntDataMap(type, _ghost); for (UInt t = 0; t < nb_tags; ++t) { std::string name(names); // std::cout << name << std::endl; names += name.size() + 1; UIntDataMap::iterator it = uint_data_map.find(name); if(it == uint_data_map.end()) { uint_data_map[name] = new Vector(0, 1, name); it = uint_data_map.find(name); } it->second->resize(nb_element); memcpy(it->second->values, data, nb_element * sizeof(UInt)); data += nb_element; it = ghost_uint_data_map.find(name); if(it == ghost_uint_data_map.end()) { ghost_uint_data_map[name] = new Vector(0, 1, name); it = ghost_uint_data_map.find(name); } it->second->resize(nb_ghost_element); memcpy(it->second->values, data, nb_ghost_element * sizeof(UInt)); data += nb_ghost_element; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildSurfaceID(Mesh & mesh) { AKANTU_DEBUG_IN(); // Vector node_offset; // Vector node_to_elem; CSR node_to_elem; /// Get list of surface elements UInt spatial_dimension = mesh.getSpatialDimension(); // buildNode2Elements(mesh, node_offset, node_to_elem, spatial_dimension-1); buildNode2Elements(mesh, node_to_elem, spatial_dimension-1); /// Find which types of elements have been linearized const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; UInt nb_types = type_list.size(); ElementType lin_element_type[nb_types]; UInt nb_lin_types = 0; UInt nb_nodes_per_element[nb_types]; UInt nb_nodes_per_element_p1[nb_types]; UInt * conn_val[nb_types]; UInt nb_element[nb_types+1]; ElementType type_p1; for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; if(Mesh::getSpatialDimension(type) != spatial_dimension) continue; ElementType facet_type = mesh.getFacetElementType(type); lin_element_type[nb_lin_types] = facet_type; nb_nodes_per_element[nb_lin_types] = Mesh::getNbNodesPerElement(facet_type); type_p1 = Mesh::getP1ElementType(facet_type); nb_nodes_per_element_p1[nb_lin_types] = Mesh::getNbNodesPerElement(type_p1); conn_val[nb_lin_types] = mesh.getConnectivity(facet_type, _not_ghost).values; nb_element[nb_lin_types] = mesh.getNbElement(facet_type, _not_ghost); nb_lin_types++; } for (UInt i = 1; i < nb_lin_types; ++i) nb_element[i] += nb_element[i+1]; for (UInt i = nb_lin_types; i > 0; --i) nb_element[i] = nb_element[i-1]; nb_element[0] = 0; /// Find close surfaces Vector surface_value_id(1, nb_element[nb_lin_types], -1); Int * surf_val = surface_value_id.values; UInt nb_surfaces = 0; UInt nb_cecked_elements; UInt nb_elements_to_ceck; UInt * elements_to_ceck = new UInt [nb_element[nb_lin_types]]; memset(elements_to_ceck, 0, nb_element[nb_lin_types]*sizeof(UInt)); for (UInt lin_el = 0; lin_el < nb_element[nb_lin_types]; ++lin_el) { if(surf_val[lin_el] != -1) continue; /* Surface id already assigned */ /* First element of new surface */ surf_val[lin_el] = nb_surfaces; nb_cecked_elements = 0; nb_elements_to_ceck = 1; memset(elements_to_ceck, 0, nb_element[nb_lin_types]*sizeof(UInt)); elements_to_ceck[0] = lin_el; // Find others elements belonging to this surface while(nb_cecked_elements < nb_elements_to_ceck) { UInt ceck_lin_el = elements_to_ceck[nb_cecked_elements]; // Transform linearized index of element into ElementType one UInt lin_type_nb = 0; while (ceck_lin_el >= nb_element[lin_type_nb+1]) lin_type_nb++; UInt ceck_el = ceck_lin_el - nb_element[lin_type_nb]; // Get connected elements UInt el_offset = ceck_el*nb_nodes_per_element[lin_type_nb]; for (UInt n = 0; n < nb_nodes_per_element_p1[lin_type_nb]; ++n) { UInt node_id = conn_val[lin_type_nb][el_offset + n]; CSR::iterator it_n; for (it_n = node_to_elem.begin(node_id); it_n != node_to_elem.end(node_id); ++it_n) { // for (UInt i = node_offset.values[node_id]; i < node_offset.values[node_id+1]; ++i) { if(surf_val[*it_n] == -1) { /* Found new surface element */ surf_val[*it_n] = nb_surfaces; elements_to_ceck[nb_elements_to_ceck] = *it_n; nb_elements_to_ceck++; } // if(surf_val[node_to_elem.values[i]] == -1) { /* Found new surface element */ // surf_val[node_to_elem.values[i]] = nb_surfaces; // elements_to_ceck[nb_elements_to_ceck] = node_to_elem.values[i]; // nb_elements_to_ceck++; // } } } nb_cecked_elements++; } nb_surfaces++; } delete [] elements_to_ceck; /// Transform local linearized element index in the global one for (UInt i = 0; i < nb_lin_types; ++i) nb_element[i] = nb_element[i+1] - nb_element[i]; UInt el_offset = 0; for (UInt type_it = 0; type_it < nb_lin_types; ++type_it) { ElementType type = lin_element_type[type_it]; Vector * surf_id_type = mesh.getSurfaceIDPointer(type, _not_ghost); surf_id_type->resize(nb_element[type_it]); surf_id_type->clear(); for (UInt el = 0; el < nb_element[type_it]; ++el) surf_id_type->values[el] = surf_val[el+el_offset]; el_offset += nb_element[type_it]; } /// Set nb_surfaces in mesh mesh.nb_surfaces = nb_surfaces; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshUtils::buildNodesPerSurface(const Mesh & mesh, CSR & nodes_per_surface) { AKANTU_DEBUG_IN(); UInt nb_surfaces = mesh.getNbSurfaces(); UInt nb_nodes = mesh.getNbNodes(); UInt spatial_dimension = mesh.getSpatialDimension(); //surface elements UInt nb_facet_types = 0; ElementType facet_type[_max_element_type]; Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator end = mesh.lastType(spatial_dimension); for(; it != end; ++it) { facet_type[nb_facet_types++] = mesh.getFacetElementType(*it); } UInt * surface_nodes_id = new UInt[nb_nodes*nb_surfaces]; std::fill_n(surface_nodes_id, nb_surfaces*nb_nodes, 0); for(UInt t = 0; t < nb_facet_types; ++t) { ElementType type = facet_type[t]; UInt nb_element = mesh.getNbElement(type); UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); UInt * connecticity = mesh.getConnectivity(type, _not_ghost).values; UInt * surface_id = mesh.getSurfaceID(type, _not_ghost).values;; for (UInt el = 0; el < nb_element; ++el) { UInt offset = *surface_id * nb_nodes; for (UInt n = 0; n < nb_nodes_per_element; ++n) { surface_nodes_id[offset + *connecticity] = 1; ++connecticity; } ++surface_id; } } nodes_per_surface.resizeRows(nb_surfaces); nodes_per_surface.clearRows(); UInt * surface_nodes_id_tmp = surface_nodes_id; for (UInt s = 0; s < nb_surfaces; ++s) for (UInt n = 0; n < nb_nodes; ++n) nodes_per_surface.rowOffset(s) += *surface_nodes_id_tmp++; nodes_per_surface.countToCSR(); nodes_per_surface.resizeCols(); nodes_per_surface.beginInsertions(); surface_nodes_id_tmp = surface_nodes_id; for (UInt s = 0; s < nb_surfaces; ++s) for (UInt n = 0; n < nb_nodes; ++n) { if (*surface_nodes_id_tmp == 1) nodes_per_surface.insertInRow(s, n); surface_nodes_id_tmp++; } nodes_per_surface.endInsertions(); delete [] surface_nodes_id; AKANTU_DEBUG_OUT(); } __END_AKANTU__ // LocalWords: ElementType diff --git a/src/model/solid_mechanics/material_inline_impl.cc b/src/model/solid_mechanics/material_inline_impl.cc index 56d04b632..b41be60ca 100644 --- a/src/model/solid_mechanics/material_inline_impl.cc +++ b/src/model/solid_mechanics/material_inline_impl.cc @@ -1,468 +1,471 @@ /** * @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 * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ __END_AKANTU__ #include "solid_mechanics_model.hh" #include __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ inline UInt Material::addElement(const ElementType & type, UInt element, const GhostType & ghost_type) { element_filter(type, ghost_type).push_back(element); return element_filter(type, ghost_type).getSize()-1; } /* -------------------------------------------------------------------------- */ inline UInt Material::getTangentStiffnessVoigtSize(UInt dim) const { return (dim * (dim - 1) / 2 + dim); } /* -------------------------------------------------------------------------- */ template inline void Material::gradUToF(const types::RMatrix & grad_u, types::RMatrix & F) { UInt size_F = F.size(); AKANTU_DEBUG_ASSERT(F.size() >= grad_u.size() && grad_u.size() == dim, "The dimension of the tensor F should be greater or equal to the dimension of the tensor grad_u."); for (UInt i = 0; i < dim; ++i) for (UInt j = 0; j < dim; ++j) F(i, j) = grad_u(i, j); for (UInt i = 0; i < size_F; ++i) F(i, i) += 1; } /* -------------------------------------------------------------------------- */ inline void Material::rightCauchy(const types::RMatrix & F, types::RMatrix & C) { C.mul(F, F); } /* -------------------------------------------------------------------------- */ inline void Material::leftCauchy(const types::RMatrix & F, types::RMatrix & B) { B.mul(F, F); } /* -------------------------------------------------------------------------- */ inline void Material::computePotentialEnergyOnQuad(types::RMatrix & grad_u, types::RMatrix & sigma, Real & epot) { epot = 0.; for (UInt i = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j) epot += sigma(i, j) * grad_u(i, j); epot *= .5; } /* -------------------------------------------------------------------------- */ template inline void Material::transferBMatrixToSymVoigtBMatrix(const types::RMatrix & B, types::RMatrix & Bvoigt, UInt nb_nodes_per_element) const { Bvoigt.clear(); for (UInt i = 0; i < dim; ++i) for (UInt n = 0; n < nb_nodes_per_element; ++n) Bvoigt(i, i + n*dim) = B(n, i); if(dim == 2) { ///in 2D, fill the @f$ [\frac{\partial N_i}{\partial x}, \frac{\partial N_i}{\partial y}]@f$ row for (UInt n = 0; n < nb_nodes_per_element; ++n) { Bvoigt(2, 1 + n*2) = B(n, 0); Bvoigt(2, 0 + n*2) = B(n, 1); } } if(dim == 3) { for (UInt n = 0; n < nb_nodes_per_element; ++n) { Real dndx = B(n, 0); Real dndy = B(n, 1); Real dndz = B(n, 2); ///in 3D, fill the @f$ [0, \frac{\partial N_i}{\partial y}, \frac{N_i}{\partial z}]@f$ row Bvoigt(3, 1 + n*3) = dndz; Bvoigt(3, 2 + n*3) = dndy; ///in 3D, fill the @f$ [\frac{\partial N_i}{\partial x}, 0, \frac{N_i}{\partial z}]@f$ row Bvoigt(4, 0 + n*3) = dndz; Bvoigt(4, 2 + n*3) = dndx; ///in 3D, fill the @f$ [\frac{\partial N_i}{\partial x}, \frac{N_i}{\partial y}, 0]@f$ row Bvoigt(5, 0 + n*3) = dndy; Bvoigt(5, 1 + n*3) = dndx; } } } /* -------------------------------------------------------------------------- */ template inline void Material::buildElementalFieldInterpolationCoodinates(__attribute__((unused)) const types::RMatrix & coordinates, __attribute__((unused)) types::RMatrix & coordMatrix) { AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ template<> inline void Material::buildElementalFieldInterpolationCoodinates<_triangle_3>(const types::RMatrix & coordinates, types::RMatrix & coordMatrix) { for (UInt i = 0; i < coordinates.rows(); ++i) coordMatrix(i, 0) = 1; } /* -------------------------------------------------------------------------- */ template<> inline void Material::buildElementalFieldInterpolationCoodinates<_triangle_6>(const types::RMatrix & coordinates, types::RMatrix & coordMatrix) { UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(_triangle_6); for (UInt i = 0; i < coordinates.rows(); ++i) { coordMatrix(i, 0) = 1; for (UInt j = 1; j < nb_quadrature_points; ++j) coordMatrix(i, j) = coordinates(i, j-1); } } /** * @todo Write a more efficient interpolation for quadrangles by * dropping unnecessary quadrature points * */ /* -------------------------------------------------------------------------- */ template<> inline void Material::buildElementalFieldInterpolationCoodinates<_quadrangle_4>(const types::RMatrix & coordinates, types::RMatrix & coordMatrix) { for (UInt i = 0; i < coordinates.rows(); ++i) { Real x = coordinates(i, 0); Real y = coordinates(i, 1); coordMatrix(i, 0) = 1; coordMatrix(i, 1) = x; coordMatrix(i, 2) = y; coordMatrix(i, 3) = x * y; } } /* -------------------------------------------------------------------------- */ template<> inline void Material::buildElementalFieldInterpolationCoodinates<_quadrangle_8>(const types::RMatrix & coordinates, types::RMatrix & coordMatrix) { for (UInt i = 0; i < coordinates.rows(); ++i) { UInt j = 0; Real x = coordinates(i, 0); Real y = coordinates(i, 1); for (UInt e = 0; e <= 2; ++e) { for (UInt n = 0; n <= 2; ++n) { coordMatrix(i, j) = std::pow(x, e) * std::pow(y, n); ++j; } } } } /* -------------------------------------------------------------------------- */ template inline UInt Material::getSizeElementalFieldInterpolationCoodinates() { return model->getFEM().getNbQuadraturePoints(type); } /* -------------------------------------------------------------------------- */ template void Material::registerParam(std::string name, T & variable, T default_value, ParamAccessType type, std::string description) { AKANTU_DEBUG_IN(); params.registerParam(name, variable, default_value, type, description); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::registerParam(std::string name, T & variable, ParamAccessType type, std::string description) { AKANTU_DEBUG_IN(); params.registerParam(name, variable, type, description); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline UInt Material::getNbDataForElements(const Vector & elements, SynchronizationTag tag) const { //UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(element.type); //UInt size = 0; if(tag == _gst_smm_stress) { return spatial_dimension * spatial_dimension * sizeof(Real) * this->getNbQuadraturePoints(elements); } return 0; } /* -------------------------------------------------------------------------- */ inline void Material::packElementData(CommunicationBuffer & buffer, const Vector & elements, SynchronizationTag tag) const { if(tag == _gst_smm_stress) { packElementDataHelper(stress, buffer, elements); // UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(element.type); // Vector::const_iterator stress_it = stress(element.type, _not_ghost).begin(spatial_dimension, spatial_dimension); // stress_it += element.element * nb_quadrature_points; // for (UInt q = 0; q < nb_quadrature_points; ++q, ++stress_it) // buffer << *stress_it; } } /* -------------------------------------------------------------------------- */ inline void Material::unpackElementData(CommunicationBuffer & buffer, const Vector & elements, SynchronizationTag tag) { if(tag == _gst_smm_stress) { unpackElementDataHelper(stress, buffer, elements); // UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(element.type); // Vector::iterator stress_it = stress(element.type, _ghost).begin(spatial_dimension, spatial_dimension); // stress_it += element.element * nb_quadrature_points; // for (UInt q = 0; q < nb_quadrature_points; ++q, ++stress_it) // buffer >> *stress_it; } } /* -------------------------------------------------------------------------- */ inline UInt Material::getNbQuadraturePoints(const Vector & elements) const { UInt nb_quad = 0; Vector::const_iterator it = elements.begin(); Vector::const_iterator end = elements.end(); for (; it != end; ++it) { const Element & el = *it; nb_quad += this->model->getFEM().getNbQuadraturePoints(el.type, el.ghost_type); } return nb_quad; } /* -------------------------------------------------------------------------- */ template inline void Material::packElementDataHelper(const ByElementTypeVector & data_to_pack, CommunicationBuffer & buffer, const Vector & elements) const { packUnpackElementDataHelper(const_cast &>(data_to_pack), buffer, elements); } /* -------------------------------------------------------------------------- */ template inline void Material::unpackElementDataHelper(ByElementTypeVector & data_to_unpack, CommunicationBuffer & buffer, const Vector & elements) const { packUnpackElementDataHelper(data_to_unpack, buffer, elements); } /* -------------------------------------------------------------------------- */ template inline void Material::packUnpackElementDataHelper(ByElementTypeVector & data_to_pack, CommunicationBuffer & buffer, const Vector & element) const { ElementType current_element_type = _not_defined; GhostType current_ghost_type = _casper; UInt nb_quad_per_elem = 0; UInt nb_component = 0; Vector * vect = NULL; Vector * element_index_material = NULL; Vector::const_iterator it = element.begin(); Vector::const_iterator end = element.end(); for (; it != end; ++it) { const Element & el = *it; if(el.type != current_element_type || el.ghost_type != current_ghost_type) { current_element_type = el.type; current_ghost_type = el.ghost_type; vect = &data_to_pack(el.type, el.ghost_type); element_index_material = &(this->model->getElementIndexByMaterial(current_element_type, current_ghost_type)); nb_quad_per_elem = this->model->getFEM().getNbQuadraturePoints(el.type, el.ghost_type); nb_component = vect->getNbComponent(); } - UInt el_id = (*element_index_material)(el.element); + UInt el_id = (*element_index_material)(el.element, 0); types::Vector data(vect->storage() + el_id * nb_component * nb_quad_per_elem, nb_component * nb_quad_per_elem); if(pack_helper) buffer << data; else buffer >> data; } } /* -------------------------------------------------------------------------- */ template inline T Material::getParam(const ID & param) const { try { return params.get(param); } catch (...) { AKANTU_EXCEPTION("No parameter " << param << " in the material " << id); } } /* -------------------------------------------------------------------------- */ template inline void Material::setParam(const ID & param, T value) { try { params.set(param, value); } catch(...) { AKANTU_EXCEPTION("No parameter " << param << " in the material " << id); } updateInternalParameters(); } /* -------------------------------------------------------------------------- */ template void Material::removeQuadraturePointsFromVectors(ByElementTypeVector & data, const ByElementTypeUInt & new_numbering) { for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; ByElementTypeVector::type_iterator it = new_numbering.firstType(0, gt, _ek_not_defined); ByElementTypeVector::type_iterator end = new_numbering.lastType(0, gt, _ek_not_defined); for (; it != end; ++it) { ElementType type = *it; if(data.exists(type, gt)){ const Vector & renumbering = new_numbering(type, gt); Vector & vect = data(type, gt); UInt nb_quad_per_elem = this->model->getFEM().getNbQuadraturePoints(type, gt); UInt nb_component = vect.getNbComponent(); Vector tmp(renumbering.getSize()*nb_quad_per_elem, nb_component); UInt new_size = 0; for (UInt i = 0; i < vect.getSize(); ++i) { UInt new_i = renumbering(i); if(new_i != UInt(-1)) { memcpy(tmp.storage() + new_i * nb_component * nb_quad_per_elem, vect.storage() + i * nb_component * nb_quad_per_elem, nb_component * nb_quad_per_elem * sizeof(T)); ++new_size; } } tmp.resize(new_size * nb_quad_per_elem); vect.copy(tmp); } } } } /* -------------------------------------------------------------------------- */ inline void Material::onElementsAdded(__attribute__((unused)) const Vector & element_list) { for (std::map::iterator it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) { resizeInternalVector(*(it->second)); } for (std::map::iterator it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) { resizeInternalVector(*(it->second)); } } /* -------------------------------------------------------------------------- */ inline void Material::onElementsRemoved(const Vector & element_list, const ByElementTypeUInt & new_numbering) { + UInt my_num = model->getInternalIndexFromID(id); + ByElementTypeUInt material_local_new_numbering("remove mat filter elem", id); Vector::const_iterator el_begin = element_list.begin(); Vector::const_iterator el_end = element_list.end(); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; ByElementTypeVector::type_iterator it = new_numbering.firstType(0, gt, _ek_not_defined); ByElementTypeVector::type_iterator end = new_numbering.lastType(0, gt, _ek_not_defined); for (; it != end; ++it) { ElementType type = *it; if(element_filter.exists(type, gt)){ Vector & elem_filter = element_filter(type, gt); Vector & element_index_material = model->getElementIndexByMaterial(type, gt); element_index_material.resize(model->getFEM().getMesh().getNbElement(type, gt)); // all materials will resize of the same size... if(!material_local_new_numbering.exists(type, gt)) material_local_new_numbering.alloc(elem_filter.getSize(), 1, type, gt); Vector & mat_renumbering = material_local_new_numbering(type, gt); const Vector & renumbering = new_numbering(type, gt); Vector elem_filter_tmp; UInt ni = 0; Element el; el.type = type; el.ghost_type = gt; for (UInt i = 0; i < elem_filter.getSize(); ++i) { el.element = elem_filter(i); if(std::find(el_begin, el_end, el) == el_end) { UInt new_el = renumbering(el.element); AKANTU_DEBUG_ASSERT(new_el != UInt(-1), "A not removed element as been badly renumbered"); elem_filter_tmp.push_back(new_el); mat_renumbering(i) = ni; - element_index_material(new_el) = ni; + element_index_material(new_el, 0) = ni; + element_index_material(new_el, 1) = my_num; ++ni; } else { mat_renumbering(i) = UInt(-1); } } elem_filter.resize(elem_filter_tmp.getSize()); elem_filter.copy(elem_filter); } } } for (std::map::iterator it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) { this->removeQuadraturePointsFromVectors(*(it->second), material_local_new_numbering); } for (std::map::iterator it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) { this->removeQuadraturePointsFromVectors(*(it->second), material_local_new_numbering); } } diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc index 42f935416..8f54dee5c 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model.cc @@ -1,1256 +1,1255 @@ /** * @file solid_mechanics_model.cc * @author Nicolas Richart * @date Thu Jul 22 14:35:38 2010 * * @brief Implementation of the SolidMechanicsModel class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_math.hh" #include "aka_common.hh" #include "solid_mechanics_model.hh" #include "integration_scheme_2nd_order.hh" #include "static_communicator.hh" #include "dof_synchronizer.hh" #ifdef AKANTU_USE_MUMPS #include "solver_mumps.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /** * A solid mechanics model need a mesh and a dimension to be created. the model * by it self can not do a lot, the good init functions should be called in * order to configure the model depending on what we want to do. * * @param mesh mesh representing the model we want to simulate * @param dim spatial dimension of the problem, if dim = 0 (default value) the * dimension of the problem is assumed to be the on of the mesh * @param id an id to identify the model */ SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : Model(id, memory_id), time_step(NAN), f_m2a(1.0), mass_matrix(NULL), velocity_damping_matrix(NULL), stiffness_matrix(NULL), jacobian_matrix(NULL), element_material("element_material", id), element_index_by_material("element index by material", id), integrator(NULL), increment_flag(false), solver(NULL), spatial_dimension(dim), mesh(mesh) { AKANTU_DEBUG_IN(); createSynchronizerRegistry(this); if (spatial_dimension == 0) spatial_dimension = mesh.getSpatialDimension(); registerFEMObject("SolidMechanicsFEM", mesh, spatial_dimension); this->displacement = NULL; this->mass = NULL; this->velocity = NULL; this->acceleration = NULL; this->force = NULL; this->residual = NULL; this->boundary = NULL; this->increment = NULL; this->increment_acceleration = NULL; this->dof_synchronizer = NULL; materials.clear(); mesh.registerEventHandler(*this); 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; if(solver) delete solver; if(mass_matrix) delete mass_matrix; if(velocity_damping_matrix) delete velocity_damping_matrix; if(stiffness_matrix && stiffness_matrix != jacobian_matrix) delete stiffness_matrix; if(jacobian_matrix) delete jacobian_matrix; if(dof_synchronizer) delete dof_synchronizer; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Initialisation */ /* -------------------------------------------------------------------------- */ /** * This function groups many of the initialization in on function. For most of * basics case the function should be enough. The functions initialize the * model, the internal vectors, set them to 0, and depending on the parameters * it also initialize the explicit or implicit solver. * * @param material_file the file containing the materials to use * @param method the analysis method wanted. See the akantu::AnalysisMethod for * the different possibilites */ void SolidMechanicsModel::initFull(std::string material_file, AnalysisMethod analysis_method) { method = analysis_method; // initialize the model initModel(); // initialize the vectors initVectors(); // set the initial condition to 0 force->clear(); velocity->clear(); acceleration->clear(); displacement->clear(); // initialize pcb if(pbc_pair.size()!=0) initPBC(); // initialize the time integration schemes switch(method) { case _explicit_dynamic: initExplicit(); break; case _implicit_dynamic: initImplicit(true); break; case _static: initImplicit(false); break; } // initialize the materials if(material_file != "") { readMaterials(material_file); initMaterials(); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initParallel(MeshPartition * partition, DataAccessor * data_accessor) { AKANTU_DEBUG_IN(); if (data_accessor == NULL) data_accessor = this; Synchronizer & synch_parallel = createParallelSynch(partition,data_accessor); synch_registry->registerSynchronizer(synch_parallel,_gst_smm_mass); // synch_registry->registerSynchronizer(synch_parallel,_gst_smm_for_strain); synch_registry->registerSynchronizer(synch_parallel,_gst_smm_stress); synch_registry->registerSynchronizer(synch_parallel,_gst_smm_boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initFEMBoundary(bool create_surface) { if(create_surface) MeshUtils::buildFacets(mesh); FEM & fem_boundary = getFEMBoundary(); fem_boundary.initShapeFunctions(); fem_boundary.computeNormalsOnControlPoints(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initExplicit() { AKANTU_DEBUG_IN(); method = _explicit_dynamic; if (integrator) delete integrator; integrator = new CentralDifference(); UInt nb_nodes = acceleration->getSize(); UInt nb_degree_of_freedom = acceleration->getNbComponent(); std::stringstream sstr; sstr << id << ":increment_acceleration"; increment_acceleration = &(alloc(sstr.str(), nb_nodes, nb_degree_of_freedom, Real())); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Allocate all the needed vectors. By default their are not necessarily set to * 0 * */ void SolidMechanicsModel::initVectors() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.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, spatial_dimension, 0)); 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"; current_position = &(alloc(sstr_curp.str(), 0, spatial_dimension, REAL_INIT_VALUE)); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_not_defined); for(; it != end; ++it) { UInt nb_element = mesh.getNbElement(*it, gt); element_material.alloc(nb_element, 1, *it, gt); } } dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian * */ void SolidMechanicsModel::initModel() { /// \todo add the current position as a parameter to initShapeFunctions for /// large deformation getFEM().initShapeFunctions(_not_ghost); getFEM().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initPBC() { Model::initPBC(); registerPBCSynchronizer(); // as long as there are ones on the diagonal of the matrix, we can put boudandary true for slaves std::map::iterator it = pbc_pair.begin(); std::map::iterator end = pbc_pair.end(); UInt dim = mesh.getSpatialDimension(); while(it != end) { for (UInt i=0; iregisterSynchronizer(*synch, _gst_smm_uv); synch_registry->registerSynchronizer(*synch, _gst_smm_mass); synch_registry->registerSynchronizer(*synch, _gst_smm_res); changeLocalEquationNumberforPBC(pbc_pair, mesh.getSpatialDimension()); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateCurrentPosition() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); current_position->resize(nb_nodes); Real * current_position_val = current_position->values; Real * position_val = mesh.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++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initializeUpdateResidualData() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); residual->resize(nb_nodes); /// copy the forces in residual for boundary conditions memcpy(residual->values, force->values, nb_nodes*spatial_dimension*sizeof(Real)); updateCurrentPosition(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Explicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /** * This function compute the second member of the motion equation. That is to * say the sum of forces @f$ r = F_{ext} - F_{int} @f$. @f$ F_{ext} @f$ is given * by the user in the force vector, and @f$ F_{int} @f$ is computed as @f$ * F_{int} = \int_{\Omega} N \sigma d\Omega@f$ * */ void SolidMechanicsModel::updateResidual(bool need_initialize) { AKANTU_DEBUG_IN(); // f = f_ext - f_int // f = f_ext if (need_initialize) initializeUpdateResidualData(); if (method == _explicit_dynamic) { // f -= fint // start synchronization synch_registry->asynchronousSynchronize(_gst_smm_uv); synch_registry->waitEndSynchronize(_gst_smm_uv); // communicate the displacement // synch_registry->asynchronousSynchronize(_gst_smm_for_strain); std::vector::iterator mat_it; // call update residual on each local elements for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllStresses(_not_ghost); } #ifdef AKANTU_DAMAGE_NON_LOCAL /* ------------------------------------------------------------------------ */ /* Computation of the non local part */ synch_registry->asynchronousSynchronize(_gst_mnl_for_average); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_not_ghost); } synch_registry->waitEndSynchronize(_gst_mnl_for_average); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_ghost); } #endif /* ------------------------------------------------------------------------ */ /* assembling the forces internal */ // communicate the strain synch_registry->asynchronousSynchronize(_gst_smm_stress); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.assembleResidual(_not_ghost); } // finalize communications synch_registry->waitEndSynchronize(_gst_smm_stress); // synch_registry->waitEndSynchronize(_gst_smm_for_strain); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.assembleResidual(_ghost); } } else { Vector * Ku = new Vector(*displacement, true, "Ku"); *Ku *= *stiffness_matrix; *residual -= *Ku; // updateSupportReaction(); delete Ku; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeStresses() { if (method == _explicit_dynamic) { // start synchronization synch_registry->asynchronousSynchronize(_gst_smm_uv); synch_registry->waitEndSynchronize(_gst_smm_uv); // compute stresses on all local elements for each materials std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllStresses(_not_ghost); } /* ------------------------------------------------------------------------ */ #ifdef AKANTU_DAMAGE_NON_LOCAL /* Computation of the non local part */ synch_registry->asynchronousSynchronize(_gst_mnl_for_average); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_not_ghost); } synch_registry->waitEndSynchronize(_gst_mnl_for_average); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_ghost); } #endif } else { std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllStressesFromTangentModuli(_not_ghost); } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateResidualInternal() { AKANTU_DEBUG_IN(); // f = f_ext - f_int - Ma - Cv = r - Ma - Cv; if(method != _static) { // f -= Ma if(mass_matrix) { // if full mass_matrix Vector * Ma = new Vector(*acceleration, true, "Ma"); *Ma *= *mass_matrix; /// \todo check unit conversion for implicit dynamics // *Ma /= f_m2a *residual -= *Ma; delete Ma; } else { // else lumped mass UInt nb_nodes = acceleration->getSize(); UInt nb_degree_of_freedom = acceleration->getNbComponent(); Real * mass_val = mass->values; Real * accel_val = acceleration->values; Real * res_val = residual->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) { if(!(*boundary_val)) { *res_val -= *accel_val * *mass_val /f_m2a; } boundary_val++; res_val++; mass_val++; accel_val++; } } // f -= Cv if(velocity_damping_matrix) { Vector * Cv = new Vector(*velocity); *Cv *= *velocity_damping_matrix; *residual -= *Cv; delete Cv; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateAcceleration() { AKANTU_DEBUG_IN(); updateResidualInternal(); if(method == _explicit_dynamic && !mass_matrix) { /* residual = residual_{n+1} - M * acceleration_n therefore solution = increment acceleration not acceleration */ solveLumped(*increment_acceleration, *mass, *residual, *boundary, f_m2a); } else { solveDynamic(*increment_acceleration); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveLumped(Vector & x, const Vector & A, const Vector & b, const Vector & boundary, Real alpha) { Real * A_val = A.storage(); Real * b_val = b.storage(); Real * x_val = x.storage(); bool * boundary_val = boundary.storage(); UInt nb_degrees_of_freedom = x.getSize() * x.getNbComponent(); for (UInt n = 0; n < nb_degrees_of_freedom; ++n) { if(!(*boundary_val)) { *x_val = alpha * (*b_val / *A_val); } x_val++; A_val++; b_val++; boundary_val++; } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitPred() { AKANTU_DEBUG_IN(); if(increment_flag) { memcpy(increment->values, displacement->values, displacement->getSize()*displacement->getNbComponent()*sizeof(Real)); } AKANTU_DEBUG_ASSERT(integrator,"itegrator should have been allocated: " << "have called initExplicit ? " << "or initImplicit ?"); integrator->integrationSchemePred(time_step, *displacement, *velocity, *acceleration, *boundary); if(increment_flag) { Real * inc_val = increment->values; Real * dis_val = displacement->values; UInt nb_nodes = displacement->getSize(); for (UInt n = 0; n < nb_nodes; ++n) { *inc_val = *dis_val - *inc_val; inc_val++; dis_val++; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitCorr() { AKANTU_DEBUG_IN(); integrator->integrationSchemeCorrAccel(time_step, *displacement, *velocity, *acceleration, *boundary, *increment_acceleration); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Implicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /** * Initialize the solver and create the sparse matrices needed. * */ void SolidMechanicsModel::initSolver(__attribute__((unused)) SolverOptions & options) { #if !defined(AKANTU_USE_MUMPS) // or other solver in the future \todo add AKANTU_HAS_SOLVER in CMake AKANTU_DEBUG_ERROR("You should at least activate one solver."); #else UInt nb_global_node = mesh.getNbGlobalNodes(); std::stringstream sstr; sstr << id << ":jacobian_matrix"; jacobian_matrix = new SparseMatrix(nb_global_node * spatial_dimension, _symmetric, spatial_dimension, sstr.str(), memory_id); // dof_synchronizer->initGlobalDOFEquationNumbers(); jacobian_matrix->buildProfile(mesh, *dof_synchronizer); #ifdef AKANTU_USE_MUMPS std::stringstream sstr_solv; sstr_solv << id << ":solver"; solver = new SolverMumps(*jacobian_matrix, sstr_solv.str()); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif //AKANTU_USE_MUMPS solver->initialize(options); #endif //AKANTU_HAS_SOLVER } /* -------------------------------------------------------------------------- */ /** * Initialize the implicit solver, either for dynamic or static cases, * * @param dynamic */ void SolidMechanicsModel::initImplicit(bool dynamic, SolverOptions & solver_options) { AKANTU_DEBUG_IN(); method = dynamic ? _implicit_dynamic : _static; initSolver(solver_options); if(method == _implicit_dynamic) { if(integrator) delete integrator; integrator = new TrapezoidalRule2(); } std::stringstream sstr; sstr << id << ":stiffness_matrix"; stiffness_matrix = new SparseMatrix(*jacobian_matrix, sstr.str(), memory_id); // } else { // stiffness_matrix = jacobian_matrix; // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initialAcceleration() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Solving Ma = f"); Solver * acc_solver = NULL; std::stringstream sstr; sstr << id << ":tmp_mass_matrix"; SparseMatrix * tmp_mass = new SparseMatrix(*mass_matrix, sstr.str(), memory_id); #ifdef AKANTU_USE_MUMPS std::stringstream sstr_solver; sstr << id << ":solver_mass_matrix"; acc_solver = new SolverMumps(*mass_matrix, sstr_solver.str()); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif //AKANTU_USE_MUMPS acc_solver->initialize(); tmp_mass->applyBoundary(*boundary); acc_solver->setRHS(*residual); acc_solver->solve(*acceleration); delete acc_solver; delete tmp_mass; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); stiffness_matrix->clear(); // call compute stiffness matrix on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->assembleStiffnessMatrix(_not_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveDynamic() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL, "You should first initialize the implicit solver and assemble the stiffness matrix"); AKANTU_DEBUG_ASSERT(mass_matrix != NULL, "You should first initialize the implicit solver and assemble the mass matrix"); if(!increment) setIncrementFlagOn(); updateResidualInternal(); solveDynamic(*increment); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveStatic() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Solving Ku = f"); AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL, "You should first initialize the implicit solver and assemble the stiffness matrix"); UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent(); // if(method != _static) jacobian_matrix->copyContent(*stiffness_matrix); jacobian_matrix->applyBoundary(*boundary); solver->setRHS(*residual); if(!increment) setIncrementFlagOn(); solver->solve(*increment); Real * increment_val = increment->values; Real * displacement_val = displacement->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) { if(!(*boundary_val)) { *displacement_val += *increment_val; } displacement_val++; boundary_val++; increment_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance) { Real error; bool tmp = testConvergenceIncrement(tolerance, error); AKANTU_DEBUG_INFO("Norm of increment : " << error); return tmp; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance, Real & error) { AKANTU_DEBUG_IN(); UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent(); error = 0; Real norm[2] = {0., 0.}; Real * increment_val = increment->storage(); bool * boundary_val = boundary->storage(); Real * displacement_val = displacement->storage(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); for (UInt d = 0; d < nb_degree_of_freedom; ++d) { if(!(*boundary_val) && is_local_node) { norm[0] += *increment_val * *increment_val; norm[1] += *displacement_val * *displacement_val; } boundary_val++; increment_val++; displacement_val++; } } StaticCommunicator::getStaticCommunicator().allReduce(norm, 2, _so_sum); norm[0] = sqrt(norm[0]); norm[1] = sqrt(norm[1]); AKANTU_DEBUG_ASSERT(!Math::isnan(norm[0]), "Something goes wrong in the solve phase"); AKANTU_DEBUG_OUT(); error = norm[0] / norm[1]; return (error < tolerance); } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceResidual(Real tolerance) { Real error; bool tmp = testConvergenceResidual(tolerance, error); AKANTU_DEBUG_INFO("Norm of residual : " << error); return tmp; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceResidual(Real tolerance, Real & norm) { AKANTU_DEBUG_IN(); UInt nb_nodes = residual->getSize(); norm = 0; Real * residual_val = residual->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); if(is_local_node) { for (UInt d = 0; d < spatial_dimension; ++d) { if(!(*boundary_val)) { norm += *residual_val * *residual_val; } boundary_val++; residual_val++; } } else { boundary_val += spatial_dimension; residual_val += spatial_dimension; } } StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum); norm = sqrt(norm); AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something goes wrong in the solve phase"); AKANTU_DEBUG_OUT(); return (norm < tolerance); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::implicitPred() { AKANTU_DEBUG_IN(); integrator->integrationSchemePred(time_step, *displacement, *velocity, *acceleration, *boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::implicitCorr() { AKANTU_DEBUG_IN(); integrator->integrationSchemeCorrDispl(time_step, *displacement, *velocity, *acceleration, *boundary, *increment); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Information */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::synchronizeBoundaries() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(synch_registry,"Synchronizer registry was not initialized." << " Did you call initParallel?"); synch_registry->synchronize(_gst_smm_boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::synchronizeResidual() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(synch_registry,"Synchronizer registry was not initialized." << " Did you call initPBC?"); synch_registry->synchronize(_gst_smm_res); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setIncrementFlagOn() { AKANTU_DEBUG_IN(); if(!increment) { UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_inc; sstr_inc << id << ":increment"; increment = &(alloc(sstr_inc.str(), nb_nodes, spatial_dimension, 0)); } increment_flag = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep() { AKANTU_DEBUG_IN(); Real min_dt = getStableTimeStep(_not_ghost); /// reduction min over all processors StaticCommunicator::getStaticCommunicator().allReduce(&min_dt, 1, _so_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); Material ** mat_val = &(materials.at(0)); Real min_dt = std::numeric_limits::max(); Real * coord = mesh.getNodes().values; Real * disp_val = displacement->values; Element elem; elem.ghost_type = ghost_type; Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator end = mesh.lastType(spatial_dimension, ghost_type); for(; it != end; ++it) { elem.type = *it; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(*it); UInt nb_element = mesh.getNbElement(*it); UInt * conn = mesh.getConnectivity(*it, ghost_type).storage(); UInt * elem_mat_val = element_material(*it, ghost_type).storage(); 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; elem.element = el; 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 = getFEM().getElementInradius(u, *it); Real el_dt = mat_val[elem_mat_val[el]]->getStableTimeStep(el_size, elem); min_dt = min_dt > el_dt ? el_dt : min_dt; } delete [] u; } AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getPotentialEnergy() { AKANTU_DEBUG_IN(); Real energy = 0.; /// call update residual on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { energy += (*mat_it)->getPotentialEnergy(); } /// reduction sum over all processors StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy() { AKANTU_DEBUG_IN(); Real ekin = 0.; UInt nb_nodes = mesh.getNbNodes(); Real * vel_val = velocity->values; Real * mass_val = mass->values; for (UInt n = 0; n < nb_nodes; ++n) { Real mv2 = 0; bool is_local_node = mesh.isLocalOrMasterNode(n); bool is_not_pbc_slave_node = !getIsPBCSlaveNode(n); bool count_node = is_local_node && is_not_pbc_slave_node; for (UInt i = 0; i < spatial_dimension; ++i) { mv2 += count_node * *vel_val * *vel_val * *mass_val; vel_val++; mass_val++; } ekin += mv2; } StaticCommunicator::getStaticCommunicator().allReduce(&ekin, 1, _so_sum); AKANTU_DEBUG_OUT(); return ekin * .5; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getExternalWork() { AKANTU_DEBUG_IN(); Real * velo = velocity->storage(); Real * forc = force->storage(); Real * resi = residual->storage(); bool * boun = boundary->storage(); Real work = 0.; UInt nb_nodes = mesh.getNbNodes(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); bool is_not_pbc_slave_node = !getIsPBCSlaveNode(n); bool count_node = is_local_node && is_not_pbc_slave_node; for (UInt i = 0; i < spatial_dimension; ++i) { if(*boun) work -= count_node * *resi * *velo * time_step; else work += count_node * *forc * *velo * time_step; ++velo; ++forc; ++resi; ++boun; } } StaticCommunicator::getStaticCommunicator().allReduce(&work, 1, _so_sum); AKANTU_DEBUG_OUT(); return work; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(std::string id) { AKANTU_DEBUG_IN(); if (id == "kinetic") { return getKineticEnergy(); } else if (id == "external work"){ return getExternalWork(); } Real energy = 0.; std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { energy += (*mat_it)->getEnergy(id); } /// reduction sum over all processors StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesAdded(const Vector & nodes_list) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); if(displacement) displacement->resize(nb_nodes); if(mass ) mass ->resize(nb_nodes); if(velocity ) velocity ->resize(nb_nodes); if(acceleration) acceleration->resize(nb_nodes); if(force ) force ->resize(nb_nodes); if(residual ) residual ->resize(nb_nodes); if(boundary ) boundary ->resize(nb_nodes); if(increment_acceleration) increment_acceleration->resize(nb_nodes); if(increment) increment->resize(nb_nodes); delete dof_synchronizer; dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onNodesAdded(nodes_list); } if(method != _explicit_dynamic) AKANTU_DEBUG_TO_IMPLEMENT(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsAdded(const Vector & element_list) { AKANTU_DEBUG_IN(); getFEM().initShapeFunctions(_not_ghost); getFEM().initShapeFunctions(_ghost); Vector::const_iterator it = element_list.begin(); Vector::const_iterator end = element_list.end(); /// \todo have rules to choose the correct material Material & mat = *materials[0]; for (; it != end; ++it) { const Element & elem = *it; element_material(elem.type, elem.ghost_type).push_back(UInt(0)); UInt mat_index = mat.addElement(elem.type, elem.element, elem.ghost_type); - element_index_by_material(elem.type, elem.ghost_type).push_back(mat_index); + UInt id[2]; + id[0] = mat_index; id[1] = 0; + element_index_by_material(elem.type, elem.ghost_type).push_back(id); } std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onElementsAdded(element_list); } if(method != _explicit_dynamic) AKANTU_DEBUG_TO_IMPLEMENT(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsRemoved(const Vector & element_list, const ByElementTypeUInt & new_numbering) { - //element_material.onElementsRemoved(new_numbering); - //element_index_by_material.onElementsRemoved(new_numbering); - std::cout << "NbNodes before purify " << mesh.getNbNodes() << std::endl; // MeshUtils::purifyMesh(mesh); std::cout << "NbNodes after purify " << mesh.getNbNodes() << std::endl; getFEM().initShapeFunctions(_not_ghost); getFEM().initShapeFunctions(_ghost); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onElementsRemoved(element_list, new_numbering); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesRemoved(const Vector & element_list, const Vector & new_numbering) { std::cout << "NbNodes in u before purify " << displacement->getSize() << " " << element_list.getSize() << " " << mesh.getNbNodes() << std::endl; if(displacement) mesh.removeNodesFromVector(*displacement, new_numbering); std::cout << "NbNodes in u after purify " << displacement->getSize() << std::endl; if(mass ) mesh.removeNodesFromVector(*mass , new_numbering); if(velocity ) mesh.removeNodesFromVector(*velocity , new_numbering); if(acceleration) mesh.removeNodesFromVector(*acceleration, new_numbering); if(force ) mesh.removeNodesFromVector(*force , new_numbering); if(residual ) mesh.removeNodesFromVector(*residual , new_numbering); if(boundary ) mesh.removeNodesFromVector(*boundary , new_numbering); if(increment_acceleration) mesh.removeNodesFromVector(*increment_acceleration, new_numbering); if(increment) mesh.removeNodesFromVector(*increment , new_numbering); delete dof_synchronizer; dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); } /* -------------------------------------------------------------------------- */ 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; getFEM().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; element_material.printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh index 421903c91..ba39032b6 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,588 +1,583 @@ /** * @file solid_mechanics_model.hh * @author Nicolas Richart * @date[creation] Thu Jul 22 11:51:06 2010 * @date[last modification] Thu Oct 14 14:00:06 2010 * * @brief Model of Solid Mechanics * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_MECHANICS_MODEL_HH__ #define __AKANTU_SOLID_MECHANICS_MODEL_HH__ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "model.hh" #include "data_accessor.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #include "integrator_cohesive.hh" #include "shape_cohesive.hh" #include "aka_types.hh" #include "integration_scheme_2nd_order.hh" #include "solver.hh" /* -------------------------------------------------------------------------- */ namespace akantu { class Material; class IntegrationScheme2ndOrder; class Contact; class SparseMatrix; } __BEGIN_AKANTU__ class SolidMechanicsModel : public Model, public DataAccessor, public MeshEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef FEMTemplate MyFEMType; typedef FEMTemplate< IntegratorCohesive, ShapeCohesive > MyFEMCohesiveType; SolidMechanicsModel(Mesh & mesh, UInt spatial_dimension = 0, const ID & id = "solid_mechanics_model", const MemoryID & memory_id = 0); virtual ~SolidMechanicsModel(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize completely the model virtual void initFull(std::string material_file, AnalysisMethod method = _explicit_dynamic); /// initialize the fem object needed for boundary conditions void initFEMBoundary(bool create_surface = true); /// register the tags associated with the parallel synchronizer void initParallel(MeshPartition * partition, DataAccessor * data_accessor=NULL); /// allocate all vectors void initVectors(); /// initialize all internal arrays for materials void initMaterials(); /// initialize the model virtual void initModel(); /// init PBC synchronizer void initPBC(); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* PBC */ /* ------------------------------------------------------------------------ */ public: /// change the equation number for proper assembly when using PBC void changeEquationNumberforPBC(std::map & pbc_pair); /// synchronize Residual for output void synchronizeResidual(); protected: /// register PBC synchronizer void registerPBCSynchronizer(); /* ------------------------------------------------------------------------ */ /* Explicit */ /* ------------------------------------------------------------------------ */ public: /// initialize the stuff for the explicit scheme void initExplicit(); /// initialize the array needed by updateResidual (residual, current_position) void initializeUpdateResidualData(); /// update the current position vector void updateCurrentPosition(); /// assemble the residual for the explicit scheme void updateResidual(bool need_initialize = true); /** * \brief compute the acceleration from the residual * this function is the explicit equivalent to solveDynamic in implicit * In the case of lumped mass just divide the residual by the mass * In the case of not lumped mass call solveDynamic<_acceleration_corrector> */ void updateAcceleration(); /// Solve the system @f[ A x = \alpha b @f] with A a lumped matrix void solveLumped(Vector & x, const Vector & A, const Vector & b, const Vector & boundary, Real alpha); /// explicit integration predictor void explicitPred(); /// explicit integration corrector void explicitCorr(); /* ------------------------------------------------------------------------ */ /* Implicit */ /* ------------------------------------------------------------------------ */ public: /// initialize the solver and the jacobian_matrix (called by initImplicit) void initSolver(SolverOptions & options = _solver_no_options); /// initialize the stuff for the implicit solver void initImplicit(bool dynamic = false, SolverOptions & solver_options = _solver_no_options); /// solve Ma = f to get the initial acceleration void initialAcceleration(); /// assemble the stiffness matrix void assembleStiffnessMatrix(); /// solve @f[ A\delta u = f_ext - f_int @f] in displacement void solveDynamic(); /// solve Ku = f void solveStatic(); /// test the convergence (norm of increment) bool testConvergenceIncrement(Real tolerance); bool testConvergenceIncrement(Real tolerance, Real & error); /// test the convergence (norm of residual) bool testConvergenceResidual(Real tolerance); bool testConvergenceResidual(Real tolerance, Real & error); /// implicit time integration predictor void implicitPred(); /// implicit time integration corrector void implicitCorr(); protected: /// finish the computation of residual to solve in increment void updateResidualInternal(); /// compute the support reaction and store it in force void updateSupportReaction(); /// compute A and solve @f[ A\delta u = f_ext - f_int @f] template void solveDynamic(Vector & increment); /* ------------------------------------------------------------------------ */ /* Explicit/Implicit */ /* ------------------------------------------------------------------------ */ public: /// compute the stresses void computeStresses(); /* ------------------------------------------------------------------------ */ /* Boundaries (solid_mechanics_model_boundary.cc) */ /* ------------------------------------------------------------------------ */ public: class SurfaceLoadFunctor { public: virtual void traction(__attribute__ ((unused)) const types::Vector & position, __attribute__ ((unused)) types::Vector & traction, __attribute__ ((unused)) const types::Vector & normal, __attribute__ ((unused)) Surface surface_id) { AKANTU_DEBUG_TO_IMPLEMENT(); } virtual void stress(__attribute__ ((unused)) const types::Vector & position, __attribute__ ((unused)) types::RMatrix & stress, __attribute__ ((unused)) const types::Vector & normal, __attribute__ ((unused)) Surface surface_id) { AKANTU_DEBUG_TO_IMPLEMENT(); } }; /// compute force vector from a function(x,y,z) that describe stresses void computeForcesFromFunction(BoundaryFunction in_function, BoundaryFunctionType function_type) __attribute__((deprecated)); template void computeForcesFromFunction(Functor & functor, BoundaryFunctionType function_type); /// integrate a force on the boundary by providing a stress tensor void computeForcesByStressTensor(const Vector & stresses, const ElementType & type, const GhostType & ghost_type); /// integrate a force on the boundary by providing a traction vector void computeForcesByTractionVector(const Vector & tractions, const ElementType & type, const GhostType & ghost_type); /// synchronize the ghost element boundaries values void synchronizeBoundaries(); /* ------------------------------------------------------------------------ */ /* Materials (solid_mechanics_model_material.cc) */ /* ------------------------------------------------------------------------ */ public: /// register a material in the dynamic database Material & registerNewMaterial(const ID & mat_type, const std::string & opt_param = ""); template Material & registerNewCustomMaterial(const ID & mat_type, const std::string & opt_param = ""); /// read the material files to instantiate all the materials void readMaterials(const std::string & filename); /// read a custom material with a keyword and class as template template UInt readCustomMaterial(const std::string & filename, const std::string & keyword); /// Use a UIntData in the mesh to specify the material to use per element void setMaterialIDsFromIntData(const std::string & data_name); -protected: - // /// read properties part of a material file and create the material - // template - // Material * readMaterialProperties(std::ifstream & infile, - // ID mat_id, - // UInt ¤t_line); - /* ------------------------------------------------------------------------ */ /* Mass (solid_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the mass matrix void assembleMass(); protected: /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumped(GhostType ghost_type); void assembleMass(GhostType ghost_type); /// fill a vector of rho void computeRho(Vector & rho, ElementType type, GhostType ghost_type); /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: inline virtual UInt getNbDataForElements(const Vector & elements, SynchronizationTag tag) const; inline virtual void packElementData(CommunicationBuffer & buffer, const Vector & elements, SynchronizationTag tag) const; inline virtual void unpackElementData(CommunicationBuffer & buffer, const Vector & elements, SynchronizationTag tag); inline virtual UInt getNbDataToPack(SynchronizationTag tag) const; inline virtual UInt getNbDataToUnpack(SynchronizationTag tag) const; inline virtual void packData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag) const; inline virtual void unpackData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag); protected: inline void splitElementByMaterial(const Vector & elements, Vector * elements_per_mat) const; template inline void packElementDataHelper(Vector & data_to_pack, CommunicationBuffer & buffer, const Vector & element) const; /* -------------------------------------------------------------------------- */ template inline void unpackElementDataHelper(Vector & data_to_unpack, CommunicationBuffer & buffer, const Vector & element) const; /* -------------------------------------------------------------------------- */ template inline void packUnpackElementDataHelper(Vector & data, CommunicationBuffer & buffer, const Vector & element) const; /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: virtual void onNodesAdded (const Vector & nodes_list); virtual void onNodesRemoved(const Vector & element_list, const Vector & new_numbering); virtual void onElementsAdded (const Vector & nodes_list); virtual void onElementsRemoved(const Vector & element_list, const ByElementTypeUInt & new_numbering); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the current value of the time step AKANTU_GET_MACRO(TimeStep, time_step, Real); /// set the value of the time step AKANTU_SET_MACRO(TimeStep, time_step, Real); /// get the value of the conversion from forces/ mass to acceleration AKANTU_GET_MACRO(F_M2A, f_m2a, Real); /// set the value of the conversion from forces/ mass to acceleration AKANTU_SET_MACRO(F_M2A, f_m2a, Real); /// get the SolidMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement, Vector &); /// get the SolidMechanicsModel::current_position vector \warn only consistent /// after a call to SolidMechanicsModel::updateCurrentPosition AKANTU_GET_MACRO(CurrentPosition, *current_position, const Vector &); /// get the SolidMechanicsModel::increment vector \warn only consistent if /// SolidMechanicsModel::setIncrementFlagOn has been called before AKANTU_GET_MACRO(Increment, *increment, Vector &); /// get the lumped SolidMechanicsModel::mass vector AKANTU_GET_MACRO(Mass, *mass, Vector &); /// get the SolidMechanicsModel::velocity vector AKANTU_GET_MACRO(Velocity, *velocity, Vector &); /// get the SolidMechanicsModel::acceleration vector, updated by /// SolidMechanicsModel::updateAcceleration AKANTU_GET_MACRO(Acceleration, *acceleration, Vector &); /// get the SolidMechanicsModel::force vector (boundary forces) AKANTU_GET_MACRO(Force, *force, Vector &); /// get the SolidMechanicsModel::residual vector, computed by /// SolidMechanicsModel::updateResidual AKANTU_GET_MACRO(Residual, *residual, Vector &); /// get the SolidMechanicsModel::boundary vector AKANTU_GET_MACRO(Boundary, *boundary, Vector &); /// get the SolidMechnicsModel::incrementAcceleration vector AKANTU_GET_MACRO(IncrementAcceleration, *increment_acceleration, Vector &); /// get the value of the SolidMechanicsModel::increment_flag AKANTU_GET_MACRO(IncrementFlag, increment_flag, bool); /// get the SolidMechanicsModel::element_material vector corresponding to a /// given akantu::ElementType AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementMaterial, element_material, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementMaterial, element_material, UInt); + AKANTU_GET_MACRO(ElementMaterial, element_material, const ByElementTypeVector &); /// vectors containing local material element index for each global element index AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementIndexByMaterial, element_index_by_material, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementIndexByMaterial, element_index_by_material, UInt); + AKANTU_GET_MACRO(ElementIndexByMaterial, element_index_by_material, const ByElementTypeVector &); /// get a particular material inline Material & getMaterial(UInt mat_index); inline const Material & getMaterial(UInt mat_index) const; /// give the number of materials inline UInt getNbMaterials() const { return materials.size(); }; /// give the material internal index from its id Int getInternalIndexFromID(const ID & id) const; /// compute the stable time step Real getStableTimeStep(); /// compute the potential energy Real getPotentialEnergy(); /// compute the kinetic energy Real getKineticEnergy(); /// compute the external work (for impose displacement, the velocity should be given too) Real getExternalWork(); /// get the energies Real getEnergy(std::string id); /// set the Contact object AKANTU_SET_MACRO(Contact, contact, Contact *); /** * @brief set the SolidMechanicsModel::increment_flag to on, the activate the * update of the SolidMechanicsModel::increment vector */ void setIncrementFlagOn(); /// get the stiffness matrix AKANTU_GET_MACRO(StiffnessMatrix, *stiffness_matrix, SparseMatrix &); /// get the mass matrix AKANTU_GET_MACRO(MassMatrix, *mass_matrix, SparseMatrix &); inline FEM & getFEMBoundary(std::string name = ""); /// get integrator AKANTU_GET_MACRO(Integrator, *integrator, const IntegrationScheme2ndOrder &); /// get access to the internal solver AKANTU_GET_MACRO(Solver, *solver, Solver &); protected: /// compute the stable time step Real getStableTimeStep(const GhostType & ghost_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// 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; /// accelerations array Vector * increment_acceleration; /// forces array Vector * force; /// residuals array Vector * residual; /// boundaries array Vector * boundary; /// array of current position used during update residual Vector * current_position; /// mass matrix SparseMatrix * mass_matrix; /// velocity damping matrix SparseMatrix * velocity_damping_matrix; /// stiffness matrix SparseMatrix * stiffness_matrix; /// jacobian matrix @f[A = cM + dD + K@f] with @f[c = \frac{1}{\beta \Delta t^2}, d = \frac{\gamma}{\beta \Delta t} @f] SparseMatrix * jacobian_matrix; /// materials of all elements ByElementTypeUInt element_material; /// vectors containing local material element index for each global element index ByElementTypeUInt element_index_by_material; /// list of used materials std::vector materials; /// integration scheme of second order used IntegrationScheme2ndOrder * integrator; /// increment of displacement Vector * increment; /// flag defining if the increment must be computed or not bool increment_flag; /// solver for implicit Solver * solver; /// object to resolve the contact Contact * contact; /// the spatial dimension UInt spatial_dimension; /// Mesh Mesh & mesh; AnalysisMethod method; }; __END_AKANTU__ /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "parser.hh" #include "material.hh" __BEGIN_AKANTU__ #include "solid_mechanics_model_tmpl.hh" #if defined (AKANTU_INCLUDE_INLINE_IMPL) # include "solid_mechanics_model_inline_impl.cc" #endif /// 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/src/model/solid_mechanics/solid_mechanics_model_inline_impl.cc b/src/model/solid_mechanics/solid_mechanics_model_inline_impl.cc index 0227f2d71..5a2cd25b2 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_inline_impl.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_inline_impl.cc @@ -1,526 +1,449 @@ /** * @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 * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ inline 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 const Material & SolidMechanicsModel::getMaterial(UInt mat_index) const { 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 FEM & SolidMechanicsModel::getFEMBoundary(std::string name) { return dynamic_cast(getFEMClassBoundary(name)); } /* -------------------------------------------------------------------------- */ inline void SolidMechanicsModel::splitElementByMaterial(const Vector & elements, Vector * elements_per_mat) const { ElementType current_element_type = _not_defined; GhostType current_ghost_type = _casper; UInt * elem_mat = NULL; Vector::const_iterator it = elements.begin(); Vector::const_iterator end = elements.end(); for (; it != end; ++it) { const Element & el = *it; if(el.type != current_element_type || el.ghost_type != current_ghost_type) { current_element_type = el.type; current_ghost_type = el.ghost_type; elem_mat = element_material(el.type, el.ghost_type).storage(); } elements_per_mat[elem_mat[el.element]].push_back(el); } } /* -------------------------------------------------------------------------- */ template inline void SolidMechanicsModel::packElementDataHelper(Vector & data_to_pack, CommunicationBuffer & buffer, const Vector & element) const { packUnpackElementDataHelper(data_to_pack, buffer, element); } /* -------------------------------------------------------------------------- */ template inline void SolidMechanicsModel::unpackElementDataHelper(Vector & data_to_unpack, CommunicationBuffer & buffer, const Vector & element) const { packUnpackElementDataHelper(data_to_unpack, buffer, element); } /* -------------------------------------------------------------------------- */ template inline void SolidMechanicsModel::packUnpackElementDataHelper(Vector & data, CommunicationBuffer & buffer, const Vector & elements) const { UInt nb_component = data.getNbComponent(); UInt nb_nodes_per_element = 0; ElementType current_element_type = _not_defined; GhostType current_ghost_type = _casper; UInt * conn = NULL; Vector::const_iterator it = elements.begin(); Vector::const_iterator end = elements.end(); for (; it != end; ++it) { const Element & el = *it; if(el.type != current_element_type || el.ghost_type != current_ghost_type) { current_element_type = el.type; current_ghost_type = el.ghost_type; conn = mesh.getConnectivity(el.type, el.ghost_type).storage(); nb_nodes_per_element = Mesh::getNbNodesPerElement(el.type); } UInt el_offset = el.element * nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt offset_conn = conn[el_offset + n]; types::Vector data_vect(data.storage() + offset_conn * nb_component, nb_component); if(pack_helper) buffer << data_vect; else buffer >> data_vect; } } } /* -------------------------------------------------------------------------- */ inline UInt SolidMechanicsModel::getNbDataForElements(const Vector & elements, SynchronizationTag tag) const { AKANTU_DEBUG_IN(); UInt size = 0; #ifndef AKANTU_NDEBUG size += elements.getSize() * spatial_dimension * sizeof(Real); /// position of the barycenter of the element (only for check) #endif UInt nb_nodes_per_element = 0; Vector::const_iterator it = elements.begin(); Vector::const_iterator end = elements.end(); for (; it != end; ++it) { const Element & el = *it; nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); } switch(tag) { case _gst_smm_mass: { size += nb_nodes_per_element * sizeof(Real) * spatial_dimension; // mass vector break; } case _gst_smm_for_strain: { size += nb_nodes_per_element * spatial_dimension * sizeof(Real); // displacement - - //UInt mat = element_material(element.type, _not_ghost)(element.element); - //size += materials[mat]->getNbDataToPack(element, tag); break; } case _gst_smm_boundary: { // force, displacement, boundary size += nb_nodes_per_element * spatial_dimension * (2 * sizeof(Real) + sizeof(bool)); break; } default: { } } Vector * elements_per_mat = new Vector[materials.size()]; this->splitElementByMaterial(elements, elements_per_mat); for (UInt i = 0; i < materials.size(); ++i) { size += materials[i]->getNbDataForElements(elements_per_mat[i], tag); } delete [] elements_per_mat; AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ inline void SolidMechanicsModel::packElementData(CommunicationBuffer & buffer, const Vector & elements, SynchronizationTag tag) const { AKANTU_DEBUG_IN(); - // GhostType ghost_type = _not_ghost; - - // UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); - // UInt el_offset = element.element * nb_nodes_per_element; - // UInt * conn = mesh.getConnectivity(element.type, ghost_type).storage(); - #ifndef AKANTU_NDEBUG Vector::const_iterator bit = elements.begin(); Vector::const_iterator bend = elements.end(); for (; bit != bend; ++bit) { const Element & element = *bit; types::RVector barycenter(spatial_dimension); mesh.getBarycenter(element.element, element.type, barycenter.storage(), element.ghost_type); buffer << barycenter; } #endif switch(tag) { case _gst_smm_mass: { packElementDataHelper(*mass, buffer, elements); - // for (UInt n = 0; n < nb_nodes_per_element; ++n) { - // UInt offset_conn = conn[el_offset + n]; - // Vector::iterator it_mass = mass->begin(spatial_dimension); - // buffer << it_mass[offset_conn]; - // } break; } case _gst_smm_for_strain: { packElementDataHelper(*displacement, buffer, elements); - // Vector::iterator it_disp = displacement->begin(spatial_dimension); - // for (UInt n = 0; n < nb_nodes_per_element; ++n) { - // UInt offset_conn = conn[el_offset + n]; - // buffer << it_disp[offset_conn]; - // } - - // UInt mat = element_material(element.type, ghost_type)(element.element); - // Element mat_element = element; - // mat_element.element = element_index_by_material(element.type, ghost_type)(element.element); - // materials[mat]->packData(buffer, mat_element, tag); break; } case _gst_smm_boundary: { packElementDataHelper(*force, buffer, elements); packElementDataHelper(*velocity, buffer, elements); packElementDataHelper(*boundary, buffer, elements); - // Vector::iterator it_force = force->begin(spatial_dimension); - // Vector::iterator it_velocity = velocity->begin(spatial_dimension); - // Vector::iterator > it_boundary = boundary->begin(spatial_dimension); - - // for (UInt n = 0; n < nb_nodes_per_element; ++n) { - // UInt offset_conn = conn[el_offset + n]; - - // buffer << it_force [offset_conn]; - // buffer << it_velocity[offset_conn]; - // buffer << it_boundary[offset_conn]; - // } break; } default: { - // UInt mat = element_material(element.type, ghost_type)(element.element); - // Element mat_element = element; - // mat_element.element = element_index_by_material(element.type, ghost_type)(element.element); - // materials[mat]->packData(buffer, mat_element, tag); - //AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } Vector * elements_per_mat = new Vector[materials.size()]; splitElementByMaterial(elements, elements_per_mat); for (UInt i = 0; i < materials.size(); ++i) { materials[i]->packElementData(buffer, elements_per_mat[i], tag); } delete [] elements_per_mat; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline void SolidMechanicsModel::unpackElementData(CommunicationBuffer & buffer, const Vector & elements, SynchronizationTag tag) { AKANTU_DEBUG_IN(); - // GhostType ghost_type = _ghost; - - // UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); - // UInt el_offset = element.element * nb_nodes_per_element; - // UInt * conn = mesh.getConnectivity(element.type, ghost_type).values; - #ifndef AKANTU_NDEBUG Vector::const_iterator bit = elements.begin(); Vector::const_iterator bend = elements.end(); for (; bit != bend; ++bit) { const Element & element = *bit; types::RVector barycenter_loc(spatial_dimension); mesh.getBarycenter(element.element, element.type, barycenter_loc.storage(), element.ghost_type); types::RVector barycenter(spatial_dimension); buffer >> barycenter; Real tolerance = 1e-15; for (UInt i = 0; i < spatial_dimension; ++i) { if(!(std::abs(barycenter(i) - barycenter_loc(i)) <= tolerance)) AKANTU_DEBUG_ERROR("Unpacking an unknown value for the element: " << element << "(barycenter[" << i << "] = " << barycenter_loc(i) << " and buffer[" << i << "] = " << barycenter(i) << ") - tag: " << tag); } } #endif switch(tag) { case _gst_smm_mass: { unpackElementDataHelper(*mass, buffer, elements); - // for (UInt n = 0; n < nb_nodes_per_element; ++n) { - // UInt offset_conn = conn[el_offset + n]; - // Vector::iterator it_mass = mass->begin(spatial_dimension); - // buffer >> it_mass[offset_conn]; - // } break; } case _gst_smm_for_strain: { unpackElementDataHelper(*displacement, buffer, elements); - // Vector::iterator it_disp = displacement->begin(spatial_dimension); - // for (UInt n = 0; n < nb_nodes_per_element; ++n) { - // UInt offset_conn = conn[el_offset + n]; - // buffer >> it_disp[offset_conn]; - // } - - // UInt mat = element_material(element.type, ghost_type)(element.element); - // Element mat_element = element; - // mat_element.element = element_index_by_material(element.type, ghost_type)(element.element); - // materials[mat]->unpackData(buffer, mat_element, tag); break; } case _gst_smm_boundary: { unpackElementDataHelper(*force, buffer, elements); unpackElementDataHelper(*velocity, buffer, elements); unpackElementDataHelper(*boundary, buffer, elements); - // Vector::iterator it_force = force->begin(spatial_dimension); - // Vector::iterator it_velocity = velocity->begin(spatial_dimension); - // Vector::iterator > it_boundary = boundary->begin(spatial_dimension); - - // for (UInt n = 0; n < nb_nodes_per_element; ++n) { - // UInt offset_conn = conn[el_offset + n]; - - // buffer >> it_force [offset_conn]; - // buffer >> it_velocity[offset_conn]; - // buffer >> it_boundary[offset_conn]; - // } break; } default: { - // UInt mat = element_material(element.type, ghost_type)(element.element); - // Element mat_element = element; - // mat_element.element = element_index_by_material(element.type, ghost_type)(element.element); - // materials[mat]->unpackData(buffer, mat_element, tag); - //AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } Vector * elements_per_mat = new Vector[materials.size()]; splitElementByMaterial(elements, elements_per_mat); for (UInt i = 0; i < materials.size(); ++i) { materials[i]->unpackElementData(buffer, elements_per_mat[i], tag); } delete [] elements_per_mat; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline UInt SolidMechanicsModel::getNbDataToPack(SynchronizationTag tag) const { AKANTU_DEBUG_IN(); UInt size = 0; // UInt nb_nodes = mesh.getNbNodes(); switch(tag) { case _gst_smm_uv: { size += sizeof(Real) * spatial_dimension * 2; break; } case _gst_smm_res: { size += sizeof(Real) * spatial_dimension; break; } case _gst_smm_mass: { size += sizeof(Real) * spatial_dimension; break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ inline UInt SolidMechanicsModel::getNbDataToUnpack(SynchronizationTag tag) const { AKANTU_DEBUG_IN(); UInt size = 0; // UInt nb_nodes = mesh.getNbNodes(); switch(tag) { case _gst_smm_uv: { size += sizeof(Real) * spatial_dimension * 2; break; } case _gst_smm_res: { size += sizeof(Real) * spatial_dimension; break; } case _gst_smm_mass: { size += sizeof(Real) * spatial_dimension; break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ inline void SolidMechanicsModel::packData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag) const { AKANTU_DEBUG_IN(); switch(tag) { case _gst_smm_uv: { Vector::iterator it_disp = displacement->begin(spatial_dimension); Vector::iterator it_velo = velocity->begin(spatial_dimension); buffer << it_disp[index]; buffer << it_velo[index]; break; } case _gst_smm_res: { Vector::iterator it_res = residual->begin(spatial_dimension); buffer << it_res[index]; break; } case _gst_smm_mass: { AKANTU_DEBUG_INFO("pack mass of node " << index << " which is " << (*mass)(index,0)); Vector::iterator it_mass = mass->begin(spatial_dimension); buffer << it_mass[index]; break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag) { AKANTU_DEBUG_IN(); switch(tag) { case _gst_smm_uv: { Vector::iterator it_disp = displacement->begin(spatial_dimension); Vector::iterator it_velo = velocity->begin(spatial_dimension); buffer >> it_disp[index]; buffer >> it_velo[index]; break; } case _gst_smm_res: { Vector::iterator it_res = residual->begin(spatial_dimension); buffer >> it_res[index]; break; } case _gst_smm_mass: { AKANTU_DEBUG_INFO("mass of node " << index << " was " << (*mass)(index,0)); Vector::iterator it_mass = mass->begin(spatial_dimension); buffer >> it_mass[index]; AKANTU_DEBUG_INFO("mass of node " << index << " is now " << (*mass)(index,0)); break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); } __END_AKANTU__ #include "sparse_matrix.hh" #include "solver.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template void SolidMechanicsModel::solveDynamic(Vector & increment) { AKANTU_DEBUG_INFO("Solving Ma + Cv + Ku = f"); NewmarkBeta * nmb_int = dynamic_cast(integrator); Real c = nmb_int->getAccelerationCoefficient(time_step); Real d = nmb_int->getVelocityCoefficient(time_step); Real e = nmb_int->getDisplacementCoefficient(time_step); // A = c M + d C + e K jacobian_matrix->clear(); if(type != NewmarkBeta::_acceleration_corrector) jacobian_matrix->add(*stiffness_matrix, e); jacobian_matrix->add(*mass_matrix, c); mass_matrix->saveMatrix("M.mtx"); if(velocity_damping_matrix) jacobian_matrix->add(*velocity_damping_matrix, d); jacobian_matrix->applyBoundary(*boundary); #ifndef AKANTU_NDEBUG if(AKANTU_DEBUG_TEST(dblDump)) jacobian_matrix->saveMatrix("J.mtx"); #endif jacobian_matrix->saveMatrix("J.mtx"); solver->setRHS(*residual); // solve A w = f solver->solve(increment); } /* -------------------------------------------------------------------------- */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_material.cc b/src/model/solid_mechanics/solid_mechanics_model_material.cc index 504c72845..6774966c1 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_material.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_material.cc @@ -1,207 +1,210 @@ /** * @file solid_mechanics_model_material.cc * @author Guillaume ANCIAUX * @date Thu Nov 25 10:48:53 2010 * * @brief * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" #include "material.hh" #include "aka_math.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #define AKANTU_INTANTIATE_MATERIAL_BY_DIM_NO_TMPL(dim, elem) \ material = \ &(registerNewCustomMaterial< BOOST_PP_ARRAY_ELEM(1, elem)< dim > >(mat_type, \ opt_param)) #define AKANTU_INTANTIATE_MATERIAL_BY_DIM_TMPL_EACH(r, data, i, elem) \ BOOST_PP_EXPR_IF(BOOST_PP_NOT_EQUAL(0, i), else ) \ if(opt_param == BOOST_PP_STRINGIZE(BOOST_PP_TUPLE_ELEM(2, 0, elem))) { \ material = \ &(registerNewCustomMaterial< BOOST_PP_ARRAY_ELEM(1, data)< BOOST_PP_ARRAY_ELEM(0, data), \ BOOST_PP_TUPLE_ELEM(2, 1, elem) > >(mat_type, opt_param)); \ } #define AKANTU_INTANTIATE_MATERIAL_BY_DIM_TMPL(dim, elem) \ BOOST_PP_SEQ_FOR_EACH_I(AKANTU_INTANTIATE_MATERIAL_BY_DIM_TMPL_EACH, \ (2, (dim, BOOST_PP_ARRAY_ELEM(1, elem))), \ BOOST_PP_ARRAY_ELEM(2, elem)) \ else { \ AKANTU_INTANTIATE_MATERIAL_BY_DIM_NO_TMPL(dim, elem); \ } #define AKANTU_INTANTIATE_MATERIAL_BY_DIM(dim, elem) \ BOOST_PP_IF(BOOST_PP_EQUAL(3, BOOST_PP_ARRAY_SIZE(elem) ), \ AKANTU_INTANTIATE_MATERIAL_BY_DIM_TMPL, \ AKANTU_INTANTIATE_MATERIAL_BY_DIM_NO_TMPL)(dim, elem) #define AKANTU_INTANTIATE_MATERIAL(elem) \ switch(spatial_dimension) { \ case 1: { AKANTU_INTANTIATE_MATERIAL_BY_DIM(1, elem); break; } \ case 2: { AKANTU_INTANTIATE_MATERIAL_BY_DIM(2, elem); break; } \ case 3: { AKANTU_INTANTIATE_MATERIAL_BY_DIM(3, elem); break; } \ } #define AKANTU_INTANTIATE_MATERIAL_IF(elem) \ if (mat_type == BOOST_PP_STRINGIZE(BOOST_PP_ARRAY_ELEM(0, elem))) { \ AKANTU_INTANTIATE_MATERIAL(elem); \ } #define AKANTU_INTANTIATE_OTHER_MATERIAL(r, data, elem) \ else AKANTU_INTANTIATE_MATERIAL_IF(elem) #define AKANTU_INSTANTIATE_MATERIALS() \ do { \ AKANTU_INTANTIATE_MATERIAL_IF(BOOST_PP_SEQ_HEAD(AKANTU_MATERIAL_LIST)) \ BOOST_PP_SEQ_FOR_EACH(AKANTU_INTANTIATE_OTHER_MATERIAL, \ _, \ BOOST_PP_SEQ_TAIL(AKANTU_MATERIAL_LIST)) \ else \ AKANTU_DEBUG_ERROR("Malformed material file : unknown material type " \ << mat_type); \ } while(0) /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::readMaterials(const std::string & filename) { Parser parser; parser.open(filename); std::string opt_param; std::string mat_type = parser.getNextSection("material", opt_param); while (mat_type != ""){ Material & mat = registerNewMaterial(mat_type, opt_param); parser.readSection(mat.getID(), mat); mat_type = parser.getNextSection("material", opt_param); } parser.close(); } /* -------------------------------------------------------------------------- */ Material & SolidMechanicsModel::registerNewMaterial(const ID & mat_type, const std::string & opt_param) { // UInt mat_count = materials.size(); // std::stringstream sstr_mat; sstr_mat << id << ":" << mat_count << ":" << mat_type; // Material * material; // ID mat_id = sstr_mat.str(); // // add all the new materials in the AKANTU_MATERIAL_LIST in the material.hh file // AKANTU_INSTANTIATE_MATERIALS(); // materials.push_back(material); Material * material = NULL; AKANTU_INSTANTIATE_MATERIALS(); return *material; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initMaterials() { AKANTU_DEBUG_ASSERT(materials.size() != 0, "No material to initialize !"); Material ** mat_val = &(materials.at(0)); /// @todo synchronize element material for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; /// fill the element filters of the materials using the element_material arrays Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt); Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt); for(; it != end; ++it) { UInt nb_element = mesh.getNbElement(*it, gt); UInt * elem_mat_val = element_material(*it, gt).storage(); - element_index_by_material.alloc(nb_element, 1, *it, gt); + element_index_by_material.alloc(nb_element, 2, *it, gt); + + Vector & el_id_by_mat = element_index_by_material(*it, gt); for (UInt el = 0; el < nb_element; ++el) { UInt index = mat_val[elem_mat_val[el]]->addElement(*it, el, gt); - element_index_by_material(*it, gt)(el) = index; + el_id_by_mat(el, 0) = index; + el_id_by_mat(el, 1) = elem_mat_val[el]; } } } std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { /// init internals properties (*mat_it)->initMaterial(); } synch_registry->synchronize(_gst_smm_init_mat); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setMaterialIDsFromIntData(const std::string & data_name) { for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt); Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt); for(; it != end; ++it) { try { const Vector & data = mesh.getUIntData(*it, data_name, gt); AKANTU_DEBUG_ASSERT(element_material.exists(*it, gt), "element_material for type (" << gt << ":" << *it << ") does not exists!"); element_material(*it, gt).copy(data); } catch(...) { AKANTU_DEBUG_ERROR("No data named " << data_name << " present in the mesh " << id << " for the element type " << *it); } } } } /* -------------------------------------------------------------------------- */ Int SolidMechanicsModel::getInternalIndexFromID(const ID & id) const { AKANTU_DEBUG_IN(); std::vector::const_iterator first = materials.begin(); std::vector::const_iterator last = materials.end(); for (; first != last; ++first) if ((*first)->getID() == id) { AKANTU_DEBUG_OUT(); return (first - materials.begin()); } AKANTU_DEBUG_OUT(); return -1; } __END_AKANTU__