diff --git a/cmake/Modules/CMakePackagesSystemPrivateFunctions.cmake b/cmake/Modules/CMakePackagesSystemPrivateFunctions.cmake index 65ab2b0ce..7b579961b 100644 --- a/cmake/Modules/CMakePackagesSystemPrivateFunctions.cmake +++ b/cmake/Modules/CMakePackagesSystemPrivateFunctions.cmake @@ -1,924 +1,924 @@ #=============================================================================== # @file CMakePackagesSystemPrivateFunctions.cmake # # @author Nicolas Richart # # @date creation: Sat Jul 18 2015 # @date last modification: Wed Jan 20 2016 # # @brief Set of macros used by the package system, internal functions # # @section LICENSE # # Copyright (©) 2015 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 . # #=============================================================================== # ============================================================================== # "Private" Accessors # ============================================================================== # ------------------------------------------------------------------------------ # Real name # ------------------------------------------------------------------------------ function(_package_get_real_name pkg_name real_name) set(${real_name} ${${pkg_name}} PARENT_SCOPE) endfunction() function(_package_set_real_name pkg_name real_name) set(${pkg_name} ${real_name} CACHE INTERNAL "" FORCE) endfunction() # ------------------------------------------------------------------------------ # Option name # ------------------------------------------------------------------------------ function(_package_declare_option pkg_name) string(TOUPPER "${PROJECT_NAME}" _project) _package_get_real_name(${pkg_name} _real_name) string(TOUPPER "${_real_name}" _u_package) _package_get_nature(${pkg_name} _nature) if(${_nature} MATCHES "internal" OR ${_nature} MATCHES "meta") set(_opt_name ${_project}_${_u_package}) elseif(${_nature} MATCHES "external") set(_opt_name ${_project}_USE_${_u_package}) else() set(_opt_name UNKNOWN_NATURE_${_project}_${_u_package}) endif() _package_set_variable(OPTION_NAME ${pkg_name} ${_opt_name}) endfunction() function(_package_get_option_name pkg_name opt_name) _package_get_variable(OPTION_NAME ${pkg_name} _opt_name) set(${opt_name} ${_opt_name} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Set if system package or compile external lib # ------------------------------------------------------------------------------ function(_package_set_system_option pkg_name default) string(TOUPPER "${PROJECT_NAME}" _project) _package_get_real_name(${pkg_name} _real_name) string(TOUPPER "${_real_name}" _u_package) option(${_project}_USE_SYSTEM_${_u_package} "Should akantu compile the third-party: \"${_real_name}\"" ${default}) mark_as_advanced(${_project}_USE_SYSTEM_${_u_package}) endfunction() function(_package_use_system pkg_name use) string(TOUPPER "${PROJECT_NAME}" _project) _package_get_real_name(${pkg_name} _real_name) string(TOUPPER "${_real_name}" _u_package) if(DEFINED ${_project}_USE_SYSTEM_${_u_package}) set(${use} ${${_project}_USE_SYSTEM_${_u_package}} PARENT_SCOPE) else() set(${use} TRUE PARENT_SCOPE) endif() endfunction() function(_package_set_system_script pkg_name script) _package_set_variable(COMPILE_SCRIPT ${pkg_name} "${script}") endfunction() function(_package_add_third_party_script_variable pkg_name var) _package_set_variable(VARIABLE_${var} ${pkg_name} "${ARGN}") set(${var} ${ARGN} PARENT_SCOPE) endfunction() function(_package_load_third_party_script pkg_name) if(${pkg_name}_COMPILE_SCRIPT) # set the stored variable get_cmake_property(_all_vars VARIABLES) foreach(_var ${_all_vars}) if(_var MATCHES "^${pkg_name}_VARIABLE_.*") string(REPLACE "${pkg_name}_VARIABLE_" "" _orig_var "${_var}") set(${_orig_var} ${${_var}}) endif() endforeach() _package_get_real_name(${pkg_name} _name) string(TOUPPER "${_name}" _u_name) _package_get_option_name(${pkg_name} _opt_name) if(${_opt_name}_VERSION) set(_version "${${_opt_name}_VERSION}") set(${_u_name}_VERSION "${_version}" CACHE INTERNAL "" FORCE) elseif(${_u_name}_VERSION) set(_version "${${_u_name}_VERSION}") endif() # load the script include(ExternalProject) include(${${pkg_name}_COMPILE_SCRIPT}) if(${_u_name}_LIBRARIES) _package_set_libraries(${pkg_name} ${${_u_name}_LIBRARIES}) list(APPEND _required_vars ${_u_name}_LIBRARIES) endif() if(${_u_name}_INCLUDE_DIR) _package_set_include_dir(${pkg_name} ${${_u_name}_INCLUDE_DIR}) list(APPEND _required_vars ${_u_name}_INCLUDE_DIR) endif() include(FindPackageHandleStandardArgs) if(CMAKE_VERSION VERSION_GREATER 2.8.12) find_package_handle_standard_args(${_name} REQUIRED_VARS ${_required_vars} VERSION_VAR _version FAIL_MESSAGE "Something was not configured by a the third-party script for ${_name}" ) else() find_package_handle_standard_args(${_name} "Something was not configured by a the third-party script for ${_name}" ${_required_vars} ) endif() endif() set(${pkg_name}_USE_SYSTEM_PREVIOUS FALSE CACHE INTERNAL "" FORCE) endfunction() # ------------------------------------------------------------------------------ # Nature # ------------------------------------------------------------------------------ function(_package_set_nature pkg_name nature) _package_set_variable(NATURE ${pkg_name} ${nature}) endfunction() function(_package_get_nature pkg_name nature) _package_get_variable(NATURE ${pkg_name} _nature "unknown") set(${nature} ${_nature} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Description # ------------------------------------------------------------------------------ function(_package_set_description pkg_name desc) _package_set_variable(DESC ${pkg_name} ${desc}) endfunction() function(_package_get_description pkg_name desc) _package_get_variable(DESC ${pkg_name} _desc "No description set for the package ${${pkg_name}} (${pkg_name})") set(${desc} ${_desc} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Package file name # ------------------------------------------------------------------------------ function(_package_set_filename pkg_name file) _package_set_variable(FILE ${pkg_name} ${file}) endfunction() function(_package_get_filename pkg_name file) _package_get_variable(FILE ${pkg_name} _file "No filename set for the package ${${pkg_name}}") set(${file} ${_file} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Source folder # ------------------------------------------------------------------------------ function(_package_set_sources_folder pkg_name src_folder) _package_set_variable(SRC_FOLDER ${pkg_name} ${src_folder}) endfunction() function(_package_get_sources_folder pkg_name src_folder) _package_get_variable(SRC_FOLDER ${pkg_name} _src_folder) set(${src_folder} ${_src_folder} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Test folder # ------------------------------------------------------------------------------ function(_package_set_tests_folder pkg_name test_folder) _package_set_variable(TEST_FOLDER ${pkg_name} ${test_folder}) endfunction() function(_package_get_tests_folder pkg_name test_folder) _package_get_variable(TEST_FOLDER ${pkg_name} _test_folder) set(${test_folder} ${_test_folder} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Manual folder # ------------------------------------------------------------------------------ function(_package_set_manual_folder pkg_name manual_folder) _package_set_variable(MANUAL_FOLDER ${pkg_name} ${manual_folder}) endfunction() function(_package_get_manual_folder pkg_name manual_folder) _package_get_variable(MANUAL_FOLDER ${pkg_name} _manual_folder) set(${manual_folder} ${_manual_folder} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Extra option for the find_package # ------------------------------------------------------------------------------ function(_package_set_find_package_extra_options pkg_name) _package_set_variable(FIND_PKG_OPTIONS ${pkg_name} ${ARGN}) endfunction() function(_package_get_find_package_extra_options pkg_name options) _package_get_variable(FIND_PKG_OPTIONS ${pkg_name} _options) set(${options} ${_options} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Compilation flags # ------------------------------------------------------------------------------ function(_package_set_compile_flags pkg_name lang) _package_set_variable(COMPILE_${lang}_FLAGS ${pkg_name} ${ARGN}) endfunction() function(_package_unset_compile_flags pkg_name lang) _package_variable_unset(COMPILE_${lang}_FLAGS ${pkg_name}) endfunction() function(_package_get_compile_flags pkg_name lang flags) _package_get_variable(COMPILE_${lang}_FLAGS ${pkg_name} _tmp_flags) string(REPLACE ";" " " _flags "${_tmp_flags}") set(${flags} "${_flags}" PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Include dir # ------------------------------------------------------------------------------ function(_package_set_include_dir pkg_name) _package_set_variable(INCLUDE_DIR ${pkg_name} ${ARGN}) endfunction() function(_package_get_include_dir pkg_name include_dir) _package_get_variable(INCLUDE_DIR ${pkg_name} _include_dir "") set(${include_dir} ${_include_dir} PARENT_SCOPE) endfunction() function(_package_add_include_dir pkg_name) _package_add_to_variable(INCLUDE_DIR ${pkg_name} ${ARGN}) endfunction() # ------------------------------------------------------------------------------ # Libraries # ------------------------------------------------------------------------------ function(_package_set_libraries pkg_name) _package_set_variable(LIBRARIES ${pkg_name} ${ARGN}) endfunction() function(_package_get_libraries pkg_name libraries) _package_get_variable(LIBRARIES ${pkg_name} _libraries "") set(${libraries} ${_libraries} PARENT_SCOPE) endfunction() function(_package_add_libraries pkg_name) _package_add_to_variable(LIBRARIES ${pkg_name} ${ARGN}) endfunction() # ------------------------------------------------------------------------------ # Extra dependencies like custom commands of ExternalProject # ------------------------------------------------------------------------------ function(_package_add_extra_dependency pkg_name) _package_add_to_variable(EXTRA_DEPENDENCY ${pkg_name} ${ARGN}) endfunction() function(_package_rm_extra_dependency pkg_name dep) _package_remove_from_variable(EXTRA_DEPENDENCY ${pkg_name} ${dep}) endfunction() function(_package_set_extra_dependencies pkg) _package_set_variable(EXTRA_DEPENDENCY ${pkg_name} ${ARGN}) endfunction() function(_package_get_extra_dependencies pkg deps) _package_get_variable(EXTRA_DEPENDENCY ${pkg_name} _deps "") set(${deps} ${_deps} PARENT_SCOPE) endfunction() function(_package_unset_extra_dependencies pkg_name) _package_variable_unset(EXTRA_DEPENDENCY ${pkg_name}) endfunction() # ------------------------------------------------------------------------------ # Activate/deactivate # ------------------------------------------------------------------------------ function(_package_activate pkg_name) _package_set_variable(STATE ${pkg_name} ON) endfunction() function(_package_deactivate pkg_name) _package_set_variable(STATE ${pkg_name} OFF) endfunction() function(_package_is_activated pkg_name act) _package_get_variable(STATE ${pkg_name} _state OFF) if(_state) set(${act} TRUE PARENT_SCOPE) else() set(${act} FALSE PARENT_SCOPE) endif() endfunction() function(_package_is_deactivated pkg_name act) _package_get_variable(STATE ${pkg_name} _state OFF) if(NOT _state) set(${act} TRUE PARENT_SCOPE) else() set(${act} FALSE PARENT_SCOPE) endif() endfunction() function(_package_unset_activated pkg_name) _package_variable_unset(STATE ${pkg_name}) endfunction() # ------------------------------------------------------------------------------ # Callbacks # ------------------------------------------------------------------------------ function(_package_on_enable_script pkg_name script) string(TOLOWER "${pkg_name}" _l_pkg_name) set(_output_file "${CMAKE_CURRENT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/${_l_pkg_name}.cmake") file(WRITE "${_output_file}" "${script}") _package_set_variable(CALLBACK_SCRIPT ${pkg_name} "${_output_file}") endfunction() function(_package_get_callback_script pkg_name filename) _package_get_variable(CALLBACK_SCRIPT ${pkg_name} _filename NOTFOUND) set(${filename} ${_filename} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Export list # ------------------------------------------------------------------------------ function(_package_add_to_export_list pkg_name) _package_add_to_variable(EXPORT_LIST ${pkg_name} ${ARGN}) endfunction() function(_package_get_export_list pkg_name export_list) _package_get_variable(EXPORT_LIST ${pkg_name} _export_list) set(${export_list} ${_export_list} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Direct dependencies # ------------------------------------------------------------------------------ function(_package_add_dependencies pkg_name) _package_add_to_variable(DEPENDENCIES ${pkg_name} ${ARGN}) endfunction() function(_package_get_dependencies pkg_name dependencies) _package_get_variable(DEPENDENCIES ${pkg_name} _dependencies) set(${dependencies} ${_dependencies} PARENT_SCOPE) endfunction() function(_package_unset_dependencies pkg_name) _package_variable_unset(DEPENDENCIES ${pkg_name}) endfunction() function(_package_remove_dependency pkg_name dep) set(_deps ${${pkg_name}_DEPENDENCIES}) _package_get_fdependencies(${dep} _fdeps) # check if this is the last reverse dependency list(LENGTH _fdeps len) list(FIND _fdeps ${pkg_name} pos) if((len EQUAL 1) AND (NOT pos EQUAL -1)) _package_get_description(${dep} _dep_desc) _package_get_option_name(${dep} _dep_option_name) set(${_dep_option_name} ${${dep}_OLD} CACHE BOOL "${_dep_desc}" FORCE) unset(${dep}_OLD CACHE) endif() # remove the pkg_name form the reverse dependency _package_remove_fdependency(${dep} ${pkg_name}) list(FIND _deps ${dep} pos) if(NOT pos EQUAL -1) list(REMOVE_AT _deps ${pos}) _package_set_variable(DEPENDENCIES ${pkg_name} ${_deps}) endif() endfunction() # ------------------------------------------------------------------------------ # Functions to handle reverse dependencies # ------------------------------------------------------------------------------ function(_package_set_rdependencies pkg_name) _package_set_variable(RDEPENDENCIES ${pkg_name} ${ARGN}) endfunction() function(_package_get_rdependencies pkg_name rdependencies) _package_get_variable(RDEPENDENCIES ${pkg_name} _rdependencies) set(${rdependencies} ${_rdependencies} PARENT_SCOPE) endfunction() function(_package_add_rdependency pkg_name rdep) # store the reverse dependency _package_add_to_variable(RDEPENDENCIES ${pkg_name} ${rdep}) endfunction() function(_package_remove_rdependency pkg_name rdep) _package_remove_from_variable(RDEPENDENCIES ${pkg_name} ${rdep}) endfunction() # ------------------------------------------------------------------------------ # Function to handle forcing dependencies (Package turn ON that enforce their # dependencies ON) # ------------------------------------------------------------------------------ function(_package_set_fdependencies pkg_name) _package_set_variable(FDEPENDENCIES ${pkg_name} ${ARGN}) endfunction() function(_package_get_fdependencies pkg_name fdependencies) _package_get_variable(FDEPENDENCIES ${pkg_name} _fdependencies) set(${fdependencies} ${_fdependencies} PARENT_SCOPE) endfunction() function(_package_add_fdependency pkg_name fdep) # store the reverse dependency _package_add_to_variable(FDEPENDENCIES ${pkg_name} ${fdep}) endfunction() function(_package_remove_fdependency pkg_name fdep) _package_remove_from_variable(FDEPENDENCIES ${pkg_name} ${fdep}) endfunction() # ------------------------------------------------------------------------------ # Exteral package system as apt rpm dependencies # ------------------------------------------------------------------------------ function(_package_set_package_system_dependency pkg system) string(TOUPPER "${_system}" _u_system) _package_set_variable(PACKAGE_SYSTEM_${_u_system} ${_pkg_name} ${ARGN}) endfunction() function(_package_get_package_system_dependency pkg system var) string(TOUPPER "${_system}" _u_system) _package_get_variable(PACKAGE_SYSTEM_${_u_system} ${_pkg_name} ${_deps}) set(${var} ${_deps} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Documentation related functions # ------------------------------------------------------------------------------ function(_package_set_documentation_files pkg_name) _package_set_variable(DOCUMENTATION_FILES ${pkg_name} ${ARGN}) endfunction() function(_package_get_documentation_files pkg_name doc_files) _package_get_variable(DOCUMENTATION_FILES ${pkg_name} _doc_files "") set(${doc_files} ${_doc_files} PARENT_SCOPE) endfunction() function(_package_set_documentation pkg_name) # \n replaced by && and \\ by ££ to avoid cache problems set(_doc_str "") foreach(_str ${ARGN}) set(_doc_str "${_doc_str}&&${_str}") endforeach() string(REPLACE "\\" "££" _doc_escaped "${_doc_str}") _package_set_variable(DOCUMENTATION ${pkg_name} "${_doc_str}") endfunction() function(_package_get_documentation pkg_name _doc) # \n replaced by && and \\ by ££ to avoid cache problems _package_get_variable(DOCUMENTATION ${pkg_name} _doc_tmp "") string(REPLACE "££" "\\" _doc_escaped "${_doc_tmp}") string(REPLACE "&&" "\n" _doc_newlines "${_doc_escaped}") set(${_doc} "${_doc_newlines}" PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Special boost thingy # ------------------------------------------------------------------------------ function(_package_set_boost_component_needed pkg_name) _package_add_to_variable(BOOST_COMPONENTS_NEEDED ${pkg_name} ${ARGN}) package_get_name(Boost _boost_pkg_name) _package_add_dependencies(${pkg_name} ${_boost_pkg_name}) endfunction() function(_package_get_boost_component_needed pkg_name needed) _package_get_variable(BOOST_COMPONENTS_NEEDED ${pkg_name} _needed) set(${needed} ${_needed} PARENT_SCOPE) endfunction() function(_package_load_boost_components) string(TOUPPER ${PROJECT_NAME} _project) _package_get_variable_for_activated(BOOST_COMPONENTS_NEEDED _boost_needed_components) package_get_name(Boost _pkg_name) if(_boost_needed_components) message(STATUS "Looking for Boost liraries: ${_boost_needed_components}") foreach(_comp ${_boost_needed_components}) find_package(Boost COMPONENTS ${_comp} QUIET) string(TOUPPER ${_comp} _u_comp) if(Boost_${_u_comp}_FOUND) message(STATUS " ${_comp}: FOUND") package_set_project_variable(BOOST_${_u_comp} TRUE) # Generate the libraries for the package _package_add_libraries(${_pkg_name} ${Boost_${_u_comp}_LIBRARY}) else() message(STATUS " ${_comp}: NOT FOUND") endif() endforeach() endif() endfunction() # ------------------------------------------------------------------------------ # get the list of source files for a given package # ------------------------------------------------------------------------------ function(_package_get_source_files pkg_name SRCS PUBLIC_HEADERS PRIVATE_HEADERS) string(TOUPPER ${PROJECT_NAME} _project) set(tmp_SRCS) set(tmp_PUBLIC_HEADERS) set(tmp_PRIVATE_HEADERS) foreach(_type SRCS PUBLIC_HEADERS PRIVATE_HEADERS) foreach(_file ${${pkg_name}_${_type}}) string(REPLACE "${CMAKE_CURRENT_SOURCE_DIR}/" "" _rel_file "${_file}") list(APPEND tmp_${_type} "${_rel_file}") endforeach() endforeach() set(${SRCS} ${tmp_SRCS} PARENT_SCOPE) set(${PUBLIC_HEADERS} ${tmp_PUBLIC_HEADERS} PARENT_SCOPE) set(${PRIVATE_HEADERS} ${tmp_PRIVATE_HEADERS} PARENT_SCOPE) endfunction() # ============================================================================== # Internal functions # ============================================================================== # ------------------------------------------------------------------------------ # Build the reverse dependencies from the dependencies # ------------------------------------------------------------------------------ function(_package_build_rdependencies) package_get_all_packages(_pkg_list) # set empty lists foreach(_pkg_name ${_pkg_list}) set(${_pkg_name}_rdeps) endforeach() # fill the dependencies list foreach(_pkg_name ${_pkg_list}) _package_get_dependencies(${_pkg_name} _deps) foreach(_dep_name ${_deps}) list(APPEND ${_dep_name}_rdeps ${_pkg_name}) endforeach() endforeach() # clean and set the reverse dependencies foreach(_pkg_name ${_pkg_list}) if(${_pkg_name}_rdeps) list(REMOVE_DUPLICATES ${_pkg_name}_rdeps) _package_set_rdependencies(${_pkg_name} ${${_pkg_name}_rdeps}) endif() endforeach() endfunction() # ------------------------------------------------------------------------------ # This function resolve the dependance order run the needed find_packages # ------------------------------------------------------------------------------ function(_package_load_packages) package_get_all_packages(_pkg_list) # Activate the dependencies of activated package and generate an ordered list # of dependencies set(ordered_loading_list) foreach(_pkg_name ${_pkg_list}) _package_load_dependencies_package(${_pkg_name} ordered_loading_list) endforeach() # Load the packages in the propoer order foreach(_pkg_name ${ordered_loading_list}) _package_get_option_name(${_pkg_name} _option_name) if(${_option_name}) _package_load_package(${_pkg_name}) else() # deactivate the packages than can already be deactivated _package_deactivate(${_pkg_name}) endif() endforeach() # generates the activated and unactivated lists of packages set(_packages_activated) set(_packages_deactivated) foreach(_pkg_name ${_pkg_list}) _package_is_activated(${_pkg_name} _active) if(_active) list(APPEND _packages_activated ${_pkg_name}) else() list(APPEND _packages_deactivated ${_pkg_name}) endif() endforeach() # generate the list usable by the calling code package_set_project_variable(ACTIVATED_PACKAGE_LIST "${_packages_activated}") package_set_project_variable(DEACTIVATED_PACKAGE_LIST "${_packages_deactivated}") endfunction() # ------------------------------------------------------------------------------ # This load an external package and recursively all its dependencies # ------------------------------------------------------------------------------ function(_package_load_dependencies_package pkg_name loading_list) # Get packages informations _package_get_option_name(${pkg_name} _pkg_option_name) _package_get_dependencies(${pkg_name} _dependencies) # handle the dependencies foreach(_dep_name ${_dependencies}) _package_get_description(${_dep_name} _dep_desc) _package_get_option_name(${_dep_name} _dep_option_name) _package_get_fdependencies(${_dep_name} _fdeps) if(${_pkg_option_name}) if("${_fdeps}" STREQUAL "") set(${_dep_name}_OLD ${${_dep_option_name}} CACHE INTERNAL "" FORCE) endif() # set the option to on set(${_dep_option_name} ON CACHE BOOL "${_dep_desc}" FORCE) # store the reverse dependency _package_add_fdependency(${_dep_name} ${pkg_name}) else() # check if this is the last reverse dependency list(LENGTH _fdeps len) list(FIND _fdeps ${pkg_name} pos) if((len EQUAL 1) AND (NOT pos EQUAL -1)) set(${_dep_option_name} ${${_dep_name}_OLD} CACHE BOOL "${_dep_desc}" FORCE) unset(${_dep_name}_OLD CACHE) endif() # remove the pkg_name form the reverse dependency _package_remove_fdependency(${_dep_name} ${pkg_name}) endif() # recusively load the dependencies _package_load_dependencies_package(${_dep_name} ${loading_list}) endforeach() # get the compile flags _package_get_compile_flags(${pkg_name} CXX _pkg_compile_flags) package_get_project_variable(NO_AUTO_COMPILE_FLAGS _no_auto_compile_flags) # if package option is on add it in the list if(${_pkg_option_name}) list(FIND ${loading_list} ${pkg_name} _pos) if(_pos EQUAL -1) set(_tmp_loading_list ${${loading_list}}) list(APPEND _tmp_loading_list ${pkg_name}) set(${loading_list} "${_tmp_loading_list}" PARENT_SCOPE) endif() #add the comilation flags if needed if(_pkg_compile_flags AND NOT _no_auto_compile_flags) add_flags(cxx ${_pkg_compile_flags}) endif() else() #remove the comilation flags if needed if(_pkg_comile_flags AND NOT _no_auto_compile_flags) remove_flags(cxx ${_pkg_compile_flags}) endif() endif() endfunction() # ------------------------------------------------------------------------------ # Load the package if it is an external one # ------------------------------------------------------------------------------ function(_package_load_package pkg_name) # load the package if it is an external _package_get_nature(${pkg_name} _nature) if(${_nature} MATCHES "external") _package_use_system(${pkg_name} _use_system) set(_activated TRUE) if(_use_system) _package_load_external_package(${pkg_name} _activated) else() _package_load_third_party_script(${pkg_name}) endif() if(_activated) _package_activate(${pkg_name}) elseif() _package_deactivate(${pkg_name}) endif() else(${_nature}) _package_activate(${pkg_name}) endif() endfunction() # ------------------------------------------------------------------------------ # Load external packages # ------------------------------------------------------------------------------ function(_package_load_external_package pkg_name activate) _package_get_real_name(${pkg_name} _real_name) string(TOUPPER ${_real_name} _u_package) if(NOT ${pkg_name}_USE_SYSTEM_PREVIOUS) #if system was off before clear the cache of preset variables get_cmake_property(_all_vars VARIABLES) foreach(_var ${_all_vars}) if(_var MATCHES "^${_u_package}_.*") unset(${_var} CACHE) endif() endforeach() set(${pkg_name}_USE_SYSTEM_PREVIOUS TRUE CACHE INTERNAL "" FORCE) endif() _package_get_find_package_extra_options(${pkg_name} _options) if(_options) cmake_parse_arguments(_opt_pkg "" "LANGUAGE" "PREFIX;FOUND;ARGS" ${_options}) if(_opt_pkg_UNPARSED_ARGUMENTS) message("You passed too many options for the find_package related to ${${pkg_name}} \"${_opt_pkg_UNPARSED_ARGUMENTS}\"") endif() endif() if(_opt_pkg_LANGUAGE) foreach(_language ${_opt_pkg_LANGUAGE}) enable_language(${_language}) endforeach() endif() # find the package find_package(${_real_name} REQUIRED ${_opt_pkg_ARGS}) # check if the package is found if(_opt_pkg_PREFIX) set(_package_prefix ${_opt_pkg_PREFIX}) else() set(_package_prefix ${_u_package}) endif() set(_act FALSE) set(_prefix_to_consider) if(DEFINED _opt_pkg_FOUND) if(${_opt_pkg_FOUND}) set(_act TRUE) set(_prefix_to_consider ${_package_prefix}) endif() else() foreach(_prefix ${_package_prefix}) if(${_prefix}_FOUND) set(_act TRUE) list(APPEND _prefix_to_consider ${_prefix}) endif() endforeach() endif() if(_act) foreach(_prefix ${_prefix_to_consider}) # Generate the include dir for the package if(DEFINED ${_prefix}_INCLUDE_DIRS) _package_set_include_dir(${pkg_name} ${${_prefix}_INCLUDE_DIRS}) elseif(DEFINED ${_prefix}_INCLUDE_DIR) _package_set_include_dir(${pkg_name} ${${_prefix}_INCLUDE_DIR}) elseif(DEFINED ${_prefix}_INCLUDE_PATH) _package_set_include_dir(${pkg_name} ${${_prefix}_INCLUDE_PATH}) endif() # Generate the libraries for the package if(DEFINED ${_prefix}_LIBRARIES) _package_set_libraries(${pkg_name} ${${_prefix}_LIBRARIES}) elseif(DEFINED ${_prefix}_LIBRARY) _package_set_libraries(${pkg_name} ${${_prefix}_LIBRARY}) endif() endforeach() _package_get_callback_script(${pkg_name} _script_file) if(_script_file) include("${_script_file}") endif() endif() set(${activate} ${_act} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Sanity check functions # ------------------------------------------------------------------------------ function(_package_check_files_exists) set(_message FALSE) package_get_all_packages(_pkg_list) foreach(_pkg_name ${_pkg_list}) set(_pkg_files ${${_pkg_name}_SRCS} ${${_pkg_name}_PUBLIC_HEADERS} ${${_pkg_name}_PRIVATE_HEADERS} ) _package_get_real_name(${_pkg_name} _real_name) foreach(_file ${_pkg_files}) if(NOT EXISTS "${_file}") if(NOT _message) set(_message TRUE) message("This file(s) is(are) present in a package but are not present on disk.") endif() message(" PACKAGE ${_real_name} FILE ${_file}") endif() endforeach() endforeach() if(_message) message(SEND_ERROR "Please check the content of your packages to correct this warnings") endif() endfunction() # ------------------------------------------------------------------------------ function(_package_check_files_registered) set(_pkg_files) package_get_all_packages(_pkg_list) # generates a file list of registered files foreach(_pkg_name ${_pkg_list}) list(APPEND _pkg_files ${${_pkg_name}_SRCS} ${${_pkg_name}_PUBLIC_HEADERS} ${${_pkg_name}_PRIVATE_HEADERS} ) endforeach() # generates the list of files in the source folders set(_all_src_files) foreach(_src_folder ${ARGN}) foreach(_ext "cc" "hh" "c" "h" "hpp") file(GLOB_RECURSE _src_files "${_src_folder}/*.${_ext}") list(APPEND _all_src_files ${_src_files}) endforeach() endforeach() if(_all_src_files) list(REMOVE_DUPLICATES _all_src_files) endif() set(_not_registerd_files) - # check only sources files ine the source folders + # check only sources files in the source folders foreach(_src_folder ${ARGN}) foreach(_file ${_all_src_files}) if("${_file}" MATCHES "${_src_folder}") list(FIND _pkg_files "${_file}" _index) if (_index EQUAL -1) list(APPEND _not_registerd_files ${_file}) endif() endif() endforeach() endforeach() if(AUTO_MOVE_UNKNOWN_FILES) file(MAKE_DIRECTORY ${PROJECT_BINARY_DIR}/unknown_files) endif() # warn the user and move the files if needed if(_not_registerd_files) if(EXISTS ${PROJECT_BINARY_DIR}/missing_files_in_packages) file(REMOVE ${PROJECT_BINARY_DIR}/missing_files_in_packages) endif() message("This files are present in the source folders but are not registered in any package") foreach(_file ${_not_registerd_files}) message(" ${_file}") if(AUTO_MOVE_UNKNOWN_FILES) get_filename_component(_file_name ${_file} NAME) file(RENAME ${_file} ${PROJECT_BINARY_DIR}/unknown_files/${_file_name}) endif() file(APPEND ${PROJECT_BINARY_DIR}/missing_files_in_packages "${_file} ") endforeach() if(AUTO_MOVE_UNKNOWN_FILES) message(SEND_ERROR "The files where moved in the followinf folder ${PROJECT_BINARY_DIR}/unknown_files\n Please register them in the good package or clean the sources") else() message(SEND_ERROR "Please register them in the good package or clean the sources") endif() endif() endfunction() # ------------------------------------------------------------------------------ diff --git a/src/model/dof_manager.cc b/src/model/dof_manager.cc index 4977f4b8b..718e452c6 100644 --- a/src/model/dof_manager.cc +++ b/src/model/dof_manager.cc @@ -1,498 +1,498 @@ /** * @file dof_manager.cc * * @author Nicolas Richart * * @date Wed Aug 12 09:52:30 2015 * * @brief Implementation of the common parts of the DOFManagers * * @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 "dof_manager.hh" #include "mesh.hh" #include "sparse_matrix.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ DOFManager::DOFManager(const ID & id, const MemoryID & memory_id) : Memory(id, memory_id), global_equation_number(0, 1, "global_equation_number"), mesh(NULL), local_system_size(0), pure_local_system_size(0), system_size(0) {} /* -------------------------------------------------------------------------- */ DOFManager::~DOFManager() { - NodesToElementsType::scalar_iterator nte_it = this->nodes_to_elements.begin(); - NodesToElementsType::scalar_iterator nte_end = this->nodes_to_elements.end(); + NodesToElements::scalar_iterator nte_it = this->nodes_to_elements.begin(); + NodesToElements::scalar_iterator nte_end = this->nodes_to_elements.end(); for (; nte_it != nte_end; ++nte_it) delete *nte_it; DOFStorage::iterator ds_it = dofs.begin(); DOFStorage::iterator ds_end = dofs.end(); for (; ds_it != ds_end; ++ds_it) delete ds_it->second; SparseMatricesMap::iterator sm_it = matrices.begin(); SparseMatricesMap::iterator sm_end = matrices.end(); for (; sm_it != sm_end; ++sm_it) delete sm_it->second; NonLinearSolversMap::iterator nls_it = non_linear_solvers.begin(); NonLinearSolversMap::iterator nls_end = non_linear_solvers.end(); for (; nls_it != nls_end; ++nls_it) delete nls_it->second; TimeStepSolversMap::iterator tss_it = time_step_solvers.begin(); TimeStepSolversMap::iterator tss_end = time_step_solvers.end(); for (; tss_it != tss_end; ++tss_it) delete tss_it->second; } /* -------------------------------------------------------------------------- */ void DOFManager::assembleElementalArrayLocalArray( const Array & elementary_vect, Array & array_assembeled, const ElementType & type, const GhostType & ghost_type, Real scale_factor, const Array & filter_elements) { AKANTU_DEBUG_IN(); UInt nb_element; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_degree_of_freedom = elementary_vect.getNbComponent() / nb_nodes_per_element; UInt * filter_it = NULL; if (filter_elements != empty_filter) { nb_element = filter_elements.getSize(); filter_it = filter_elements.storage(); } else { nb_element = this->mesh->getNbElement(type, ghost_type); } AKANTU_DEBUG_ASSERT(elementary_vect.getSize() == nb_element, "The vector elementary_vect(" << elementary_vect.getID() << ") has not the good size."); const Array connectivity = this->mesh->getConnectivity(type, ghost_type); Array::const_vector_iterator conn_begin = connectivity.begin(nb_nodes_per_element); Array::const_vector_iterator conn_it = conn_begin; UInt size_mat = nb_nodes_per_element * nb_degree_of_freedom; Array::const_vector_iterator elem_it = elementary_vect.begin(size_mat); for (UInt el = 0; el < nb_element; ++el, ++elem_it) { if (filter_it != NULL) conn_it = conn_begin + *filter_it; for (UInt n = 0, ld = 0; n < nb_nodes_per_element; ++n) { UInt offset_node = (*conn_it)(n)*nb_degree_of_freedom; for (UInt d = 0; d < nb_degree_of_freedom; ++d, ++ld) { array_assembeled[offset_node + d] += scale_factor * (*elem_it)(ld); } } if (filter_it != NULL) ++filter_it; else ++conn_it; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -void DOFManager::assembleElementalArrayResidual( +void DOFManager::assembleElementalArrayToResidual( const ID & dof_id, const Array & elementary_vect, const ElementType & type, const GhostType & ghost_type, Real scale_factor, const Array & filter_elements) { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_degree_of_freedom = elementary_vect.getNbComponent() / nb_nodes_per_element; Array array_localy_assembeled(this->mesh->getNbNodes(), nb_degree_of_freedom); + array_localy_assembeled.clear(); + this->assembleElementalArrayLocalArray( elementary_vect, array_localy_assembeled, type, ghost_type, scale_factor, filter_elements); this->assembleToResidual(dof_id, array_localy_assembeled, scale_factor); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::registerMesh(Mesh & mesh) { this->mesh = &mesh; this->mesh->registerEventHandler(*this, 20); UInt nb_nodes = this->mesh->getNbNodes(); this->nodes_to_elements.resize(nb_nodes); for (UInt n = 0; n < nb_nodes; ++n) { this->nodes_to_elements[n] = new std::set(); } } /* -------------------------------------------------------------------------- */ -DOFManager::DOFData::DOFData() : solution(0, 1, "solution") { } +DOFManager::DOFData::DOFData() : solution(0, 1, "solution") {} /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFs(const ID & dof_id, Array & dofs_array, const DOFSupportType & support_type) { DOFStorage::iterator it = this->dofs.find(dof_id); if (it != this->dofs.end()) { AKANTU_EXCEPTION("This dof array has already been registered"); } DOFData * dofs_storage = new DOFData(); dofs_storage->dof = &dofs_array; dofs_storage->blocked_dofs = NULL; dofs_storage->support_type = support_type; UInt nb_local_dofs = 0; UInt nb_pure_local = 0; switch (support_type) { case _dst_nodal: { nb_local_dofs = this->mesh->getNbNodes(); AKANTU_DEBUG_ASSERT( dofs_array.getSize() == nb_local_dofs, "The array of dof is too shot to be associated to nodes."); UInt nb_nodes = this->mesh->getNbNodes(); for (UInt n = 0; n < nb_nodes; ++n) { nb_pure_local += this->mesh->isLocalOrMasterNode(n) ? 1 : 0; } nb_pure_local *= dofs_array.getNbComponent(); nb_local_dofs *= dofs_array.getNbComponent(); break; } case _dst_generic: { nb_local_dofs = nb_pure_local = dofs_array.getSize() * dofs_array.getNbComponent(); break; } default: { AKANTU_EXCEPTION("This type of dofs is not handled yet."); } } this->pure_local_system_size += nb_pure_local; this->local_system_size += nb_local_dofs; StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); comm.allReduce(nb_pure_local, _so_sum); this->system_size += nb_pure_local; this->dofs[dof_id] = dofs_storage; } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFsDerivative(const ID & dof_id, UInt order, Array & dofs_derivative) { DOFStorage::iterator it = this->dofs.find(dof_id); if (it == this->dofs.end()) { AKANTU_EXCEPTION("The dof array corresponding to this derivatives has not " << "been registered yet"); } DOFData & dof = *it->second; std::vector *> & derivatives = dof.dof_derivatives; if (derivatives.size() < order) { derivatives.resize(order, NULL); } else { if (derivatives[order - 1] != NULL) { AKANTU_EXCEPTION("The dof derivatives of order " << order << " already been registered for this dof (" << dof_id << ")"); } } derivatives[order - 1] = &dofs_derivative; } /* -------------------------------------------------------------------------- */ void DOFManager::registerBlockedDOFs(const ID & dof_id, Array & blocked_dofs) { DOFStorage::iterator it = this->dofs.find(dof_id); if (it == this->dofs.end()) { AKANTU_EXCEPTION("The dof array corresponding to this derivatives has not " << "been registered yet"); } DOFData & dof = *it->second; if (dof.blocked_dofs != NULL) { AKANTU_EXCEPTION("The blocked dofs array for " << dof_id << " has already been registered"); } dof.blocked_dofs = &blocked_dofs; } /* -------------------------------------------------------------------------- */ -void DOFManager::clearJacobian() { - this->getMatrix("J").clear(); -} +void DOFManager::clearJacobian() { this->getMatrix("J").clear(); } /* -------------------------------------------------------------------------- */ void DOFManager::splitSolutionPerDOFs() { DOFStorage::iterator it = this->dofs.begin(); DOFStorage::iterator end = this->dofs.end(); for (; it != end; ++it) { DOFData & dof_data = *it->second; dof_data.solution.resize(dof_data.dof->getSize() * dof_data.dof->getNbComponent()); this->getSolutionPerDOFs(it->first, dof_data.solution); } } /* -------------------------------------------------------------------------- */ void DOFManager::registerSparseMatrix(const ID & matrix_id, SparseMatrix & matrix) { SparseMatricesMap::const_iterator it = this->matrices.find(matrix_id); if (it != this->matrices.end()) { AKANTU_EXCEPTION("The matrix " << matrix_id << " already exists in " << this->id); } this->matrices[matrix_id] = &matrix; } /* -------------------------------------------------------------------------- */ /// Get an instance of a new SparseMatrix Array & DOFManager::getNewLumpedMatrix(const ID & id) { ID matrix_id = this->id + ":lumpmtx:" + id; LumpedMatricesMap::const_iterator it = this->lumped_matrices.find(matrix_id); if (it != this->lumped_matrices.end()) { AKANTU_EXCEPTION("The lumped matrix " << matrix_id << " already exists in " << this->id); } Array & mtx = this->alloc(matrix_id, this->local_system_size, 1); this->lumped_matrices[matrix_id] = &mtx; return mtx; } /* -------------------------------------------------------------------------- */ void DOFManager::registerNonLinearSolver(const ID & non_linear_solver_id, NonLinearSolver & non_linear_solver) { NonLinearSolversMap::const_iterator it = this->non_linear_solvers.find(non_linear_solver_id); if (it != this->non_linear_solvers.end()) { AKANTU_EXCEPTION("The non linear solver " << non_linear_solver_id << " already exists in " << this->id); } this->non_linear_solvers[non_linear_solver_id] = &non_linear_solver; } /* -------------------------------------------------------------------------- */ void DOFManager::registerTimeStepSolver(const ID & time_step_solver_id, TimeStepSolver & time_step_solver) { TimeStepSolversMap::const_iterator it = this->time_step_solvers.find(time_step_solver_id); if (it != this->time_step_solvers.end()) { AKANTU_EXCEPTION("The non linear solver " << time_step_solver_id << " already exists in " << this->id); } this->time_step_solvers[time_step_solver_id] = &time_step_solver; } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManager::getMatrix(const ID & id) { ID matrix_id = this->id + ":mtx:" + id; SparseMatricesMap::const_iterator it = this->matrices.find(matrix_id); if (it == this->matrices.end()) { - AKANTU_EXCEPTION("The matrix " << matrix_id << " does not exists in " + AKANTU_SILENT_EXCEPTION("The matrix " << matrix_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ const Array & DOFManager::getLumpedMatrix(const ID & id) { ID matrix_id = this->id + ":lumpmtx:" + id; LumpedMatricesMap::const_iterator it = this->lumped_matrices.find(matrix_id); if (it == this->lumped_matrices.end()) { - AKANTU_EXCEPTION("The lumped matrix " << matrix_id << " does not exists in " + AKANTU_SILENT_EXCEPTION("The lumped matrix " << matrix_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ NonLinearSolver & DOFManager::getNonLinearSolver(const ID & id) { ID non_linear_solver_id = this->id + ":nls:" + id; NonLinearSolversMap::const_iterator it = this->non_linear_solvers.find(non_linear_solver_id); if (it == this->non_linear_solvers.end()) { AKANTU_EXCEPTION("The non linear solver " << non_linear_solver_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ TimeStepSolver & DOFManager::getTimeStepSolver(const ID & id) { ID time_step_solver_id = this->id + ":tss:" + id; TimeStepSolversMap::const_iterator it = this->time_step_solvers.find(time_step_solver_id); if (it == this->time_step_solvers.end()) { AKANTU_EXCEPTION("The non linear solver " << time_step_solver_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ void DOFManager::fillNodesToElements() { UInt spatial_dimension = this->mesh->getSpatialDimension(); Element e; UInt nb_nodes = this->mesh->getNbNodes(); for (UInt n = 0; n < nb_nodes; ++n) { this->nodes_to_elements[n]->clear(); } for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { Mesh::type_iterator first = this->mesh->firstType(spatial_dimension, *gt, _ek_not_defined); Mesh::type_iterator last = this->mesh->lastType(spatial_dimension, *gt, _ek_not_defined); e.ghost_type = *gt; for (; first != last; ++first) { ElementType type = *first; e.type = type; e.kind = Mesh::getKind(type); UInt nb_element = this->mesh->getNbElement(type, *gt); Array::const_iterator > conn_it = this->mesh->getConnectivity(type, *gt) .begin(Mesh::getNbNodesPerElement(type)); for (UInt el = 0; el < nb_element; ++el, ++conn_it) { e.element = el; const Vector & conn = *conn_it; for (UInt n = 0; n < conn.size(); ++n) nodes_to_elements[conn(n)]->insert(e); } } } } /* -------------------------------------------------------------------------- */ /* Mesh Events */ /* -------------------------------------------------------------------------- */ void DOFManager::onNodesAdded(const Array & nodes_list, __attribute__((unused)) const NewNodesEvent & event) { Array::const_scalar_iterator it = nodes_list.begin(); Array::const_scalar_iterator end = nodes_list.end(); UInt nb_nodes = this->mesh->getNbNodes(); this->nodes_to_elements.resize(nb_nodes); for (; it != end; ++it) { this->nodes_to_elements[*it] = new std::set(); } } /* -------------------------------------------------------------------------- */ void DOFManager::onNodesRemoved(const Array & nodes_list, const Array & new_numbering, __attribute__((unused)) const RemovedNodesEvent & event) { Array::const_scalar_iterator it = nodes_list.begin(); Array::const_scalar_iterator end = nodes_list.end(); for (; it != end; ++it) { delete this->nodes_to_elements[*it]; } this->mesh->removeNodesFromArray(this->nodes_to_elements, new_numbering); } /* -------------------------------------------------------------------------- */ void DOFManager::onElementsAdded(const Array & elements_list, __attribute__((unused)) const NewElementsEvent & event) { Array::const_scalar_iterator it = elements_list.begin(); Array::const_scalar_iterator end = elements_list.end(); for (; it != end; ++it) { const Element & elem = *it; if (this->mesh->getSpatialDimension(elem.type) != this->mesh->getSpatialDimension()) continue; const Array & conn = this->mesh->getConnectivity(elem.type, elem.ghost_type); UInt nb_nodes_per_elem = this->mesh->getNbNodesPerElement(elem.type); for (UInt n = 0; n < nb_nodes_per_elem; ++n) { nodes_to_elements[conn(elem.element, n)]->insert(elem); } } } /* -------------------------------------------------------------------------- */ void DOFManager::onElementsRemoved( __attribute__((unused)) const Array & elements_list, __attribute__((unused)) const ElementTypeMapArray & new_numbering, __attribute__((unused)) const RemovedElementsEvent & event) { this->fillNodesToElements(); } /* -------------------------------------------------------------------------- */ void DOFManager::onElementsChanged( __attribute__((unused)) const Array & old_elements_list, __attribute__((unused)) const Array & new_elements_list, __attribute__((unused)) const ElementTypeMapArray & new_numbering, __attribute__((unused)) const ChangedElementsEvent & event) { this->fillNodesToElements(); } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/dof_manager.hh b/src/model/dof_manager.hh index fdda5e4a3..d2af3ae6a 100644 --- a/src/model/dof_manager.hh +++ b/src/model/dof_manager.hh @@ -1,355 +1,366 @@ /** * @file dof_manager.hh * * @author Nicolas Richart * * @date Wed Jul 22 11:43:43 2015 * * @brief Class handling the different types of dofs * * @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_memory.hh" #include "non_linear_solver.hh" #include "time_step_solver.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_DOF_MANAGER_HH__ #define __AKANTU_DOF_MANAGER_HH__ __BEGIN_AKANTU__ class DOFManager : protected Memory, protected MeshEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - DOFManager(const ID & id = "dof_manager", - const MemoryID & memory_id = 0); + DOFManager(const ID & id = "dof_manager", const MemoryID & memory_id = 0); virtual ~DOFManager(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// register a mesh for dof that have a support type on nodes virtual void registerMesh(Mesh & mesh); /// register an array of degree of freedom virtual void registerDOFs(const ID & dof_id, Array & dofs_array, const DOFSupportType & support_type); /// register an array of derivatives for a particular dof array virtual void registerDOFsDerivative(const ID & dof_id, UInt order, Array & dofs_derivative); /// register array representing the blocked degree of freedoms - virtual void registerBlockedDOFs(const ID & dof_id, Array & blocked_dofs); + virtual void registerBlockedDOFs(const ID & dof_id, + Array & blocked_dofs); /// Assemble an array to the global residual array virtual void assembleToResidual(const ID & dof_id, const Array & array_to_assemble, Real scale_factor = 1.) = 0; /** * Assemble elementary values to a local array of the size nb_nodes * * nb_dof_per_node. The dof number is implicitly considered as * conn(el, n) * nb_nodes_per_element + d. * With 0 < n < nb_nodes_per_element and 0 < d < nb_dof_per_node **/ virtual void assembleElementalArrayLocalArray( const Array & elementary_vect, Array & array_assembeled, const ElementType & type, const GhostType & ghost_type, Real scale_factor = 1., const Array & filter_elements = empty_filter); /** * Assemble elementary values to the global residual array. The dof number is * implicitly considered as conn(el, n) * nb_nodes_per_element + d. * With 0 < n < nb_nodes_per_element and 0 < d < nb_dof_per_node **/ - virtual void assembleElementalArrayResidual( + virtual void assembleElementalArrayToResidual( const ID & dof_id, const Array & elementary_vect, const ElementType & type, const GhostType & ghost_type, Real scale_factor = 1., const Array & filter_elements = empty_filter); /** * Assemble elementary values to the global residual array. The dof number is * implicitly considered as conn(el, n) * nb_nodes_per_element + d. With 0 < * n < nb_nodes_per_element and 0 < d < nb_dof_per_node **/ virtual void assembleElementalMatricesToMatrix( const ID & matrix_id, const ID & dof_id, const Array & elementary_mat, const ElementType & type, - const GhostType & ghost_type, + const GhostType & ghost_type = _not_ghost, const MatrixType & elemental_matrix_type = _symmetric, const Array & filter_elements = empty_filter) = 0; + /// multiply a vector by a matrix and assemble the result to the residual + virtual void assembleMatMulVectToResidual(const ID & dof_id, const ID & A_id, + const Array x, + Real scale_factor = 1) = 0; + + /// multiply a vector by a lumped matrix and assemble the result to the + /// residual + virtual void assembleLumpedMatMulVectToResidual(const ID & dof_id, + const ID & A_id, + const Array x, + Real scale_factor = 1) = 0; + /// notation fully defined yet... // virtual void assemblePreassembledMatrix(const ID & matrix_id, // const ID & dof_id_m, // const ID & dof_id_n, // const Matrix & matrix) = 0; /// sets the residual to 0 virtual void clearResidual() = 0; /// sets the jacobian matrix to 0 - virtual void clearJacobian(); + virtual void clearJacobian() = 0; /// splits the solution storage from a global view to the per dof storages void splitSolutionPerDOFs(); protected: /// minimum functionality to implement per derived version of the DOFManager /// to allow the splitSolutionPerDOFs function to work virtual void getSolutionPerDOFs(const ID & dof_id, Array & solution_array) = 0; /// fill a Vector with the equation numbers corresponding to the given /// connectivity inline void extractElementEquationNumber( const Array & equation_numbers, const Vector & connectivity, UInt nb_degree_of_freedom, Vector & local_equation_number); /// converts local equation numbers to global equation numbers; template inline void localToGlobalEquationNumber(S & inout); /* ------------------------------------------------------------------------ */ /// register a matrix void registerSparseMatrix(const ID & matrix_id, SparseMatrix & matrix); /// register a non linear solver instantiated by a derived class void registerNonLinearSolver(const ID & non_linear_solver_id, NonLinearSolver & non_linear_solver); /// register a time step solver instantiated by a derived class void registerTimeStepSolver(const ID & time_step_solver_id, TimeStepSolver & time_step_solver); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ protected: struct DOFData; inline DOFData & getDOFData(const ID & dof_id); inline const DOFData & getDOFData(const ID & dof_id) const; public: /// get the equation numbers (in local numbering) corresponding to a dof ID inline const Array & getLocalEquationNumbers(const ID & dof_id) const; /// return the local index of the global equation number inline UInt globalToLocalEquationNumber(UInt global) const; /// Global number of dofs AKANTU_GET_MACRO(SystemSize, this->system_size, UInt); /// Local number of dofs AKANTU_GET_MACRO(LocalSystemSize, this->local_system_size, UInt); /* ------------------------------------------------------------------------ */ /* DOFs and derivatives accessors */ /* ------------------------------------------------------------------------ */ /// Get a reference to the registered dof array for a given id inline Array & getDOFs(const ID & dofs_id); /// Get a reference to the registered dof derivatives array for a given id inline Array & getDOFsDerivatives(const ID & dofs_id, UInt order); /// Get a reference to the blocked dofs array registered for the given id inline const Array & getBlockedDOFs(const ID & dofs_id) const; /// Get a reference to the solution array registered for the given id inline const Array & getSolution(const ID & dofs_id) const; /* ------------------------------------------------------------------------ */ /* Matrices accessors */ /* ------------------------------------------------------------------------ */ /// Get an instance of a new SparseMatrix virtual SparseMatrix & getNewMatrix(const ID & matrix_id, const MatrixType & matrix_type) = 0; /// Get an instance of a new SparseMatrix as a copy of the SparseMatrix /// matrix_to_copy_id virtual SparseMatrix & getNewMatrix(const ID & matrix_id, const ID & matrix_to_copy_id) = 0; /// Get the reference of an existing matrix SparseMatrix & getMatrix(const ID & matrix_id); /// Get an instance of a new lumped matrix virtual Array & getNewLumpedMatrix(const ID & matrix_id); /// Get the lumped version of a given matrix const Array & getLumpedMatrix(const ID & matrix_id); /* ------------------------------------------------------------------------ */ /* Non linear system solver */ /* ------------------------------------------------------------------------ */ /// Get instance of a non linear solver virtual NonLinearSolver & getNewNonLinearSolver( const ID & nls_solver_id, const NonLinearSolverType & _non_linear_solver_type) = 0; /// get instance of a non linear solver virtual NonLinearSolver & getNonLinearSolver(const ID & nls_solver_id); /* ------------------------------------------------------------------------ */ /* Time-Step Solver */ /* ------------------------------------------------------------------------ */ /// Get instance of a time step solver virtual TimeStepSolver & getNewTimeStepSolver(const ID & time_step_solver_id, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver) = 0; /// get instance of a time step solver virtual TimeStepSolver & getTimeStepSolver(const ID & time_step_solver_id); /* ------------------------------------------------------------------------ */ /* MeshEventHandler interface */ /* ------------------------------------------------------------------------ */ private: /// fills the nodes_to_elements structure void fillNodesToElements(); public: /// function to implement to react on akantu::NewNodesEvent virtual void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event); /// function to implement to react on akantu::RemovedNodesEvent virtual void onNodesRemoved(const Array & nodes_list, const Array & new_numbering, const RemovedNodesEvent & event); /// function to implement to react on akantu::NewElementsEvent virtual void onElementsAdded(const Array & elements_list, const NewElementsEvent & event); /// function to implement to react on akantu::RemovedElementsEvent virtual void onElementsRemoved(const Array & elements_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event); /// function to implement to react on akantu::ChangedElementsEvent virtual void onElementsChanged(const Array & old_elements_list, const Array & new_elements_list, const ElementTypeMapArray & new_numbering, const ChangedElementsEvent & event); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// dof representations in the dof manager struct DOFData { DOFData(); /// DOF support type (nodal, general) this is needed to determine how the /// dof are shared among processors DOFSupportType support_type; /// Degree of freedom array Array * dof; /// Blocked degree of freedoms array Array * blocked_dofs; /// Solution associated to the dof Array solution; /* ---------------------------------------------------------------------- */ /* data for dynamic simulations */ /* ---------------------------------------------------------------------- */ /// Degree of freedom derivatives arrays std::vector *> dof_derivatives; /// local numbering equation numbers Array local_equation_number; }; - - typedef Array< std::set * > NodesToElementsType; + typedef Array *> NodesToElements; /// This info is stored to simplify the dynamic changes - NodesToElementsType nodes_to_elements; + NodesToElements nodes_to_elements; /// equation number in global numbering Array global_equation_number; typedef unordered_map::type equation_numbers_map; /// dual information of global_equation_number equation_numbers_map global_to_local_mapping; /// type to store dofs information typedef std::map DOFStorage; /// type to store all the matrices typedef std::map SparseMatricesMap; /// type to store all the lumped matrices typedef std::map *> LumpedMatricesMap; /// type to store all the non linear solver typedef std::map NonLinearSolversMap; /// type to store all the time step solver typedef std::map TimeStepSolversMap; /// store a reference to the dof arrays DOFStorage dofs; /// list of sparse matrices that where created SparseMatricesMap matrices; /// list of lumped matrices LumpedMatricesMap lumped_matrices; /// non linear solvers storage NonLinearSolversMap non_linear_solvers; /// time step solvers storage TimeStepSolversMap time_step_solvers; /// reference to the underlying mesh Mesh * mesh; /// Total number of degrees of freedom (size with the ghosts) UInt local_system_size; /// Number of purely local dofs (size without the ghosts) UInt pure_local_system_size; /// Total number of degrees of freedom UInt system_size; }; __END_AKANTU__ #include "dof_manager_inline_impl.cc" #endif /* __AKANTU_DOF_MANAGER_HH__ */ diff --git a/src/model/dof_manager_default.cc b/src/model/dof_manager_default.cc index 582f077d7..5ff15321f 100644 --- a/src/model/dof_manager_default.cc +++ b/src/model/dof_manager_default.cc @@ -1,396 +1,505 @@ /** * @file dof_manager_default.cc * * @author Nicolas Richart * * @date Tue Aug 11 16:21:01 2015 * * @brief Implementation of the default DOFManager * * @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 "dof_manager_default.hh" #include "sparse_matrix_aij.hh" #include "time_step_solver_default.hh" #include "static_communicator.hh" #include "non_linear_solver_default.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ inline void DOFManagerDefault::addSymmetricElementalMatrixToSymmetric( SparseMatrixAIJ & matrix, const Matrix & elementary_mat, const Vector & equation_numbers, UInt max_size) { for (UInt i = 0; i < elementary_mat.rows(); ++i) { UInt c_irn = equation_numbers(i); if (c_irn < max_size) { for (UInt j = i; j < elementary_mat.cols(); ++j) { UInt c_jcn = equation_numbers(j); if (c_jcn < max_size) { matrix(c_irn, c_jcn) += elementary_mat(i, j); } } } } } /* -------------------------------------------------------------------------- */ inline void DOFManagerDefault::addUnsymmetricElementalMatrixToSymmetric( SparseMatrixAIJ & matrix, const Matrix & elementary_mat, const Vector & equation_numbers, UInt max_size) { for (UInt i = 0; i < elementary_mat.rows(); ++i) { UInt c_irn = equation_numbers(i); if (c_irn < max_size) { for (UInt j = 0; j < elementary_mat.cols(); ++j) { UInt c_jcn = equation_numbers(j); if (c_jcn < max_size) { if (c_jcn >= c_irn) { matrix(c_irn, c_jcn) += elementary_mat(i, j); } } } } } } /* -------------------------------------------------------------------------- */ inline void DOFManagerDefault::addElementalMatrixToUnsymmetric( SparseMatrixAIJ & matrix, const Matrix & elementary_mat, const Vector & equation_numbers, UInt max_size) { for (UInt i = 0; i < elementary_mat.rows(); ++i) { UInt c_irn = equation_numbers(i); if (c_irn < max_size) { for (UInt j = 0; j < elementary_mat.cols(); ++j) { UInt c_jcn = equation_numbers(j); if (c_jcn < max_size) { matrix(c_irn, c_jcn) += elementary_mat(i, j); } } } } } /* -------------------------------------------------------------------------- */ DOFManagerDefault::DOFManagerDefault(const ID & id, const MemoryID & memory_id) : DOFManager(id, memory_id), residual(0, 1, std::string(id + ":residual")), global_solution(0, 1, std::string(id + ":global_solution")), global_blocked_dofs(0, 1, std::string(id + ":global_blocked_dofs")), dofs_type(0, 1, std::string(id + ":dofs_type")) {} /* -------------------------------------------------------------------------- */ DOFManagerDefault::~DOFManagerDefault() {} /* -------------------------------------------------------------------------- */ template void DOFManagerDefault::assembleToGlobalArray( const ID & dof_id, const Array & array_to_assemble, Array & global_array, T scale_factor) { AKANTU_DEBUG_IN(); const Array & equation_number = this->getLocalEquationNumbers(dof_id); UInt nb_degree_of_freedoms = array_to_assemble.getSize() * array_to_assemble.getNbComponent(); AKANTU_DEBUG_ASSERT(equation_number.getSize() == nb_degree_of_freedoms, "The array to assemble does not have a correct size." << " (" << array_to_assemble.getID() << ")"); typename Array::const_scalar_iterator arr_it = array_to_assemble.begin_reinterpret(nb_degree_of_freedoms); Array::const_scalar_iterator equ_it = equation_number.begin(); for (UInt d = 0; d < nb_degree_of_freedoms; ++d, ++arr_it, ++equ_it) { global_array(*equ_it) += scale_factor * (*arr_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::registerDOFs(const ID & dof_id, Array & dofs_array, const DOFSupportType & support_type) { // stores the current numbers of dofs UInt local_nb_dofs = this->local_system_size; UInt pure_local_nb_dofs = this->pure_local_system_size; // update or create the dof_data DOFManager::registerDOFs(dof_id, dofs_array, support_type); // Count the number of pure local dofs per proc StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); UInt prank = comm.whoAmI(); UInt psize = comm.getNbProc(); Array nb_dofs_per_proc(psize); nb_dofs_per_proc(prank) = this->pure_local_system_size - pure_local_nb_dofs; comm.allGather(nb_dofs_per_proc); UInt first_global_dofs_id = std::accumulate( nb_dofs_per_proc.begin(), nb_dofs_per_proc.begin() + prank, 0); // nb local dofs to account for UInt nb_dofs = this->local_system_size - local_nb_dofs; DOFData & dof_data = *dofs[dof_id]; this->global_equation_number.resize(this->local_system_size); // set the equation numbers UInt first_dof_id = local_nb_dofs; dof_data.local_equation_number.resize(nb_dofs); this->dofs_type.resize(local_system_size); for (UInt d = 0; d < nb_dofs; ++d) { UInt local_eq_num = first_dof_id + d; dof_data.local_equation_number(d) = local_eq_num; UInt global_eq_num = first_global_dofs_id + d; this->global_equation_number(local_eq_num) = global_eq_num; this->global_to_local_mapping[global_eq_num] = local_eq_num; switch (support_type) { case _dst_nodal: { UInt node = d / dof_data.dof->getNbComponent(); this->dofs_type(local_eq_num) = this->mesh->getNodeType(node); break; } case _dst_generic: { this->dofs_type(local_eq_num) = _nt_normal; break; } default: { AKANTU_EXCEPTION("This type of dofs is not handled yet."); } } } this->residual.resize(this->local_system_size); this->global_solution.resize(this->local_system_size); this->global_blocked_dofs.resize(this->local_system_size); } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManagerDefault::getNewMatrix(const ID & id, const MatrixType & matrix_type) { ID matrix_id = this->id + ":mtx:" + id; SparseMatrix * sm = new SparseMatrixAIJ(*this, matrix_type, matrix_id); this->registerSparseMatrix(matrix_id, *sm); return *sm; } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManagerDefault::getNewMatrix(const ID & id, const ID & matrix_to_copy_id) { ID matrix_id = this->id + ":mtx:" + id; SparseMatrixAIJ & sm_to_copy = this->getMatrix(matrix_to_copy_id); SparseMatrix * sm = new SparseMatrixAIJ(sm_to_copy, matrix_id); this->registerSparseMatrix(matrix_id, *sm); return *sm; } /* -------------------------------------------------------------------------- */ SparseMatrixAIJ & DOFManagerDefault::getMatrix(const ID & id) { SparseMatrix & matrix = DOFManager::getMatrix(id); return dynamic_cast(matrix); } /* -------------------------------------------------------------------------- */ NonLinearSolver & DOFManagerDefault::getNewNonLinearSolver(const ID & id, const NonLinearSolverType & type) { ID non_linear_solver_id = this->id + ":nls:" + id; NonLinearSolver * nls = NULL; switch (type) { case _nls_newton_raphson: case _nls_newton_raphson_modified: { nls = new NonLinearSolverNewtonRaphson(*this, type, non_linear_solver_id, this->memory_id); break; } case _nls_linear: { nls = new NonLinearSolverLinear(*this, type, non_linear_solver_id, this->memory_id); break; } case _nls_lumped: { nls = new NonLinearSolverLumped(*this, type, non_linear_solver_id, this->memory_id); break; } default: AKANTU_EXCEPTION("The asked type of non linear solver is not supported by " "this dof manager"); } this->registerNonLinearSolver(non_linear_solver_id, *nls); return *nls; } /* -------------------------------------------------------------------------- */ TimeStepSolver & DOFManagerDefault::getNewTimeStepSolver(const ID & id, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver) { ID time_step_solver_id = this->id + ":tss:" + id; TimeStepSolver * tss = new TimeStepSolverDefault( *this, type, non_linear_solver, time_step_solver_id, this->memory_id); this->registerTimeStepSolver(time_step_solver_id, *tss); return *tss; } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::clearResidual() { this->residual.resize(this->local_system_size); this->residual.clear(); } +/* -------------------------------------------------------------------------- */ +void DOFManagerDefault::clearJacobian() { + this->getMatrix("J").clear(); +} + /* -------------------------------------------------------------------------- */ void DOFManagerDefault::updateGlobalBlockedDofs() { DOFStorage::iterator it = this->dofs.begin(); DOFStorage::iterator end = this->dofs.end(); this->global_blocked_dofs.resize(this->local_system_size); this->global_blocked_dofs.clear(); for (; it != end; ++it) { DOFData & dof_data = *it->second; this->assembleToGlobalArray(it->first, *dof_data.blocked_dofs, this->global_blocked_dofs, true); } } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::getSolutionPerDOFs(const ID & dof_id, Array & solution_array) { AKANTU_DEBUG_IN(); const Array & equation_number = this->getLocalEquationNumbers(dof_id); UInt nb_degree_of_freedoms = equation_number.getSize(); solution_array.resize(nb_degree_of_freedoms); Real * sol_it = solution_array.storage(); UInt * equ_it = equation_number.storage(); for (UInt d = 0; d < nb_degree_of_freedoms; ++d, ++sol_it, ++equ_it) { (*sol_it) = this->global_solution(*equ_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleToResidual( const ID & dof_id, const Array & array_to_assemble, Real scale_factor) { AKANTU_DEBUG_IN(); this->assembleToGlobalArray(dof_id, array_to_assemble, this->residual, scale_factor); AKANTU_DEBUG_OUT(); } +/* -------------------------------------------------------------------------- */ +void DOFManagerDefault::assembleMatMulVectToResidual(const ID & dof_id, + const ID & A_id, + const Array x, + Real scale_factor) { + SparseMatrixAIJ & A = this->getMatrix(A_id); + Array global_x(this->local_system_size, 1, 0.); + + this->assembleToGlobalArray(dof_id, x, global_x, 1.); + + A.matVecMul(global_x, this->residual, scale_factor, 1.); +} + +/* -------------------------------------------------------------------------- */ +void DOFManagerDefault::assembleLumpedMatMulVectToResidual(const ID & dof_id, + const ID & A_id, + const Array x, + Real scale_factor) { + const Array & A = this->getLumpedMatrix(A_id); + + Array global_x(this->local_system_size, 1, 0.); + this->assembleToGlobalArray(dof_id, x, global_x, scale_factor); + + Array::const_scalar_iterator A_it = A.begin(); + Array::const_scalar_iterator A_end = A.end(); + Array::const_scalar_iterator x_it = global_x.begin(); + Array::scalar_iterator r_it = this->residual.begin(); + + for (; A_it != A_end; ++A_it, ++x_it, ++r_it) { + *r_it += *A_it * *x_it; + } +} + /* -------------------------------------------------------------------------- */ void DOFManagerDefault::assembleElementalMatricesToMatrix( const ID & matrix_id, const ID & dof_id, const Array & elementary_mat, const ElementType & type, const GhostType & ghost_type, const MatrixType & elemental_matrix_type, const Array & filter_elements) { AKANTU_DEBUG_IN(); + this->addToProfile(matrix_id, dof_id); + const Array & equation_number = this->getLocalEquationNumbers(dof_id); SparseMatrixAIJ & A = this->getMatrix(matrix_id); UInt nb_element; if (ghost_type == _not_ghost) { nb_element = this->mesh->getNbElement(type); } else { AKANTU_DEBUG_TO_IMPLEMENT(); } UInt * filter_it = NULL; if (filter_elements != empty_filter) { nb_element = filter_elements.getSize(); filter_it = filter_elements.storage(); } else { nb_element = this->mesh->getNbElement(type, ghost_type); } AKANTU_DEBUG_ASSERT(elementary_mat.getSize() == nb_element, "The vector elementary_mat(" << elementary_mat.getID() << ") has not the good size."); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); - UInt nb_degree_of_freedom = elementary_mat.getNbComponent() / - (nb_nodes_per_element * nb_nodes_per_element); + + UInt nb_degree_of_freedom = this->getDOFs(dof_id).getNbComponent(); + // UInt nb_degree_of_freedom = elementary_mat.getNbComponent() / + // (nb_nodes_per_element * nb_nodes_per_element); const Array connectivity = this->mesh->getConnectivity(type, ghost_type); Array::const_vector_iterator conn_begin = connectivity.begin(nb_nodes_per_element); Array::const_vector_iterator conn_it = conn_begin; UInt size_mat = nb_nodes_per_element * nb_degree_of_freedom; Vector element_eq_nb(nb_degree_of_freedom * nb_nodes_per_element); Array::const_matrix_iterator el_mat_it = elementary_mat.begin(size_mat, size_mat); for (UInt e = 0; e < nb_element; ++e, ++el_mat_it) { if (filter_it != NULL) conn_it = conn_begin + *filter_it; this->extractElementEquationNumber(equation_number, *conn_it, nb_degree_of_freedom, element_eq_nb); this->localToGlobalEquationNumber(element_eq_nb); if (filter_it != NULL) ++filter_it; else ++conn_it; if (A.getMatrixType() == _symmetric) if (elemental_matrix_type == _symmetric) this->addSymmetricElementalMatrixToSymmetric( A, *el_mat_it, element_eq_nb, A.getSize()); else this->addUnsymmetricElementalMatrixToSymmetric( A, *el_mat_it, element_eq_nb, A.getSize()); else this->addElementalMatrixToUnsymmetric(A, *el_mat_it, element_eq_nb, A.getSize()); } AKANTU_DEBUG_OUT(); } +/* -------------------------------------------------------------------------- */ +void DOFManagerDefault::addToProfile(const ID & matrix_id, const ID & dof_id) { + AKANTU_DEBUG_IN(); + + const DOFData & dof_data = this->getDOFData(dof_id); + + if (dof_data.support_type != _dst_nodal) + return; + + std::pair mat_dof = std::make_pair(matrix_id, dof_id); + if (this->matrix_profiled_dofs.find(mat_dof) != + this->matrix_profiled_dofs.end()) + return; + + UInt nb_degree_of_freedom_per_node = dof_data.dof->getNbComponent(); + + const Array & equation_number = this->getLocalEquationNumbers(dof_id); + + SparseMatrixAIJ & A = this->getMatrix(matrix_id); + + // if(irn_jcn_to_k) delete irn_jcn_to_k; + // irn_jcn_to_k = new std::map, UInt>; + // A.clearProfile(); + + UInt size = A.getSize(); + + Mesh::type_iterator it = this->mesh->firstType( + this->mesh->getSpatialDimension(), _not_ghost, _ek_not_defined); + Mesh::type_iterator end = this->mesh->lastType( + this->mesh->getSpatialDimension(), _not_ghost, _ek_not_defined); + for (; it != end; ++it) { + UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); + + const Array & connectivity = + this->mesh->getConnectivity(*it, _not_ghost); + Array::const_vector_iterator cit = + connectivity.begin(nb_nodes_per_element); + Array::const_vector_iterator cend = + connectivity.end(nb_nodes_per_element); + + UInt size_mat = nb_nodes_per_element * nb_degree_of_freedom_per_node; + Vector element_eq_nb(size_mat); + + for (; cit != cend; ++cit) { + this->extractElementEquationNumber( + equation_number, *cit, nb_degree_of_freedom_per_node, element_eq_nb); + this->localToGlobalEquationNumber(element_eq_nb); + + for (UInt i = 0; i < size_mat; ++i) { + UInt c_irn = element_eq_nb(i); + if (c_irn < size) { + for (UInt j = 0; j < size_mat; ++j) { + UInt c_jcn = element_eq_nb(j); + if (c_jcn < size) { + A.addToProfile(c_irn, c_jcn); + } + } + } + } + } + } + + this->matrix_profiled_dofs.insert(mat_dof); + + AKANTU_DEBUG_OUT(); +} + /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/dof_manager_default.hh b/src/model/dof_manager_default.hh index b3f3cdf93..6c858889b 100644 --- a/src/model/dof_manager_default.hh +++ b/src/model/dof_manager_default.hh @@ -1,191 +1,211 @@ /** * @file dof_manager_default.hh * * @author Nicolas Richart * * @date Tue Aug 11 14:06:18 2015 * * @brief Default implementation of the dof manager * * @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 "dof_manager.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_DOF_MANAGER_DEFAULT_HH__ #define __AKANTU_DOF_MANAGER_DEFAULT_HH__ namespace akantu { class SparseMatrixAIJ; class NonLinearSolverDefault; class TimeStepSolverDefault; } __BEGIN_AKANTU__ class DOFManagerDefault : public DOFManager { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: DOFManagerDefault(const ID & id = "dof_manager_default", const MemoryID & memory_id = 0); virtual ~DOFManagerDefault(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// register an array of degree of freedom virtual void registerDOFs(const ID & dof_id, Array & dofs_array, const DOFSupportType & support_type); /// Assemble an array to the global residual array virtual void assembleToResidual(const ID & dof_id, const Array & array_to_assemble, Real scale_factor = 1.); /** * Assemble elementary values to the global matrix. The dof number is * implicitly considered as conn(el, n) * nb_nodes_per_element + d. * With 0 < n < nb_nodes_per_element and 0 < d < nb_dof_per_node **/ virtual void assembleElementalMatricesToMatrix( const ID & matrix_id, const ID & dof_id, const Array & elementary_mat, const ElementType & type, const GhostType & ghost_type, const MatrixType & elemental_matrix_type, const Array & filter_elements); + /// multiply a vector by a matrix and assemble the result to the residual + virtual void assembleMatMulVectToResidual(const ID & dof_id, const ID & A_id, + const Array x, + Real scale_factor = 1); + + /// multiply a vector by a lumped matrix and assemble the result to the + /// residual + virtual void assembleLumpedMatMulVectToResidual(const ID & dof_id, + const ID & A_id, + const Array x, + Real scale_factor = 1); + protected: /// Assemble an array to the global residual array - template + template void assembleToGlobalArray(const ID & dof_id, const Array & array_to_assemble, - Array & global_array, - T scale_factor); + Array & global_array, T scale_factor); public: /// clear the residual virtual void clearResidual(); + /// clear the jacobian matrix + virtual void clearJacobian(); + /// update the global dofs vector virtual void updateGlobalBlockedDofs(); protected: /// Get the part of the solution corresponding to the dof_id virtual void getSolutionPerDOFs(const ID & dof_id, Array & solution_array); private: /// Add a symmetric matrices to a symmetric sparse matrix void addSymmetricElementalMatrixToSymmetric( SparseMatrixAIJ & matrix, const Matrix & element_mat, const Vector & equation_numbers, UInt max_size); /// Add a unsymmetric matrices to a symmetric sparse matrix (i.e. cohesive /// elements) void addUnsymmetricElementalMatrixToSymmetric( SparseMatrixAIJ & matrix, const Matrix & element_mat, const Vector & equation_numbers, UInt max_size); /// Add a matrices to a unsymmetric sparse matrix void addElementalMatrixToUnsymmetric(SparseMatrixAIJ & matrix, const Matrix & element_mat, const Vector & equation_numbers, UInt max_size); + void addToProfile(const ID & matrix_id, const ID & dof_id); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// Get an instance of a new SparseMatrix virtual SparseMatrix & getNewMatrix(const ID & matrix_id, const MatrixType & matrix_type); /// Get an instance of a new SparseMatrix as a copy of the SparseMatrix /// matrix_to_copy_id virtual SparseMatrix & getNewMatrix(const ID & matrix_id, const ID & matrix_to_copy_id); /// Get the reference of an existing matrix SparseMatrixAIJ & getMatrix(const ID & matrix_id); /* ------------------------------------------------------------------------ */ /* Non Linear Solver */ /* ------------------------------------------------------------------------ */ /// Get instance of a non linear solver virtual NonLinearSolver & getNewNonLinearSolver(const ID & nls_solver_id, const NonLinearSolverType & _non_linear_solver_type); /* ------------------------------------------------------------------------ */ /* Time-Step Solver */ /* ------------------------------------------------------------------------ */ /// Get instance of a time step solver TimeStepSolver & getNewTimeStepSolver(const ID & id, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver); /* ------------------------------------------------------------------------ */ /// Get the solution array AKANTU_GET_MACRO_NOT_CONST(GlobalSolution, global_solution, Array &); /// Get the residual array AKANTU_GET_MACRO_NOT_CONST(Residual, residual, Array &); /// Get the blocked dofs array AKANTU_GET_MACRO(GlobalBlockedDOFs, global_blocked_dofs, const Array &); /// Get the location type of a given dof inline bool isLocalOrMasterDOF(UInt dof_num); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: typedef std::map AIJMatrixMap; typedef std::map DefaultNonLinearSolversMap; typedef std::map DefaultTimeStepSolversMap; + typedef std::set > DOFToMatrixProfile; + + /// contains the the dofs that where added to the profile of a given matrix. + DOFToMatrixProfile matrix_profiled_dofs; + /// rhs to the system of equation corresponding to the residual linked to the /// different dofs Array residual; /// solution of the system of equation corresponding to the different dofs Array global_solution; /// blocked degree of freedom in the system equation corresponding to the /// different dofs Array global_blocked_dofs; /// define the dofs type, local, shared, ghost Array dofs_type; /// Map of the different matrices stored in the dof manager AIJMatrixMap aij_matrices; /// Map of the different time step solvers stored with there real type DefaultTimeStepSolversMap default_time_step_solver_map; }; __END_AKANTU__ #include "dof_manager_default_inline_impl.cc" #endif /* __AKANTU_DOF_MANAGER_DEFAULT_HH__ */ diff --git a/src/model/dof_manager_inline_impl.cc b/src/model/dof_manager_inline_impl.cc index 8ad675b85..de43ec4e1 100644 --- a/src/model/dof_manager_inline_impl.cc +++ b/src/model/dof_manager_inline_impl.cc @@ -1,131 +1,131 @@ /** * @file dof_manager_inline_impl.cc * * @author Nicolas Richart * * @date Wed Aug 12 11:07:01 2015 * * @brief * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_DOF_MANAGER_INLINE_IMPL_CC__ #define __AKANTU_DOF_MANAGER_INLINE_IMPL_CC__ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ inline void DOFManager::extractElementEquationNumber( const Array & equation_numbers, const Vector & connectivity, UInt nb_degree_of_freedom, Vector & element_equation_number) { for (UInt i = 0, ld = 0; i < connectivity.size(); ++i) { UInt n = connectivity(i); for (UInt d = 0; d < nb_degree_of_freedom; ++d, ++ld) { element_equation_number(ld) = equation_numbers(n * nb_degree_of_freedom + d); } } } /* -------------------------------------------------------------------------- */ template inline void DOFManager::localToGlobalEquationNumber(S & inout) { for (UInt i = 0; i < inout.size(); ++i) { inout(i) = this->global_equation_number(inout(i)); } } template<> inline void DOFManager::localToGlobalEquationNumber(UInt & inout) { inout = this->global_equation_number(inout); } /* -------------------------------------------------------------------------- */ inline UInt DOFManager::globalToLocalEquationNumber(UInt global) const { equation_numbers_map::const_iterator it = this->global_to_local_mapping.find(global); AKANTU_DEBUG_ASSERT(it != this->global_to_local_mapping.end(), "This global equation number " << global << " does not exists in " << this->id); return it->second; } /* -------------------------------------------------------------------------- */ inline DOFManager::DOFData & DOFManager::getDOFData(const ID & dof_id) { DOFStorage::iterator it = this->dofs.find(dof_id); if (it == this->dofs.end()) { AKANTU_EXCEPTION("The dof " << dof_id << " does not exists in " << this->id); } return *it->second; } /* -------------------------------------------------------------------------- */ const DOFManager::DOFData & DOFManager::getDOFData(const ID & dof_id) const { DOFStorage::const_iterator it = this->dofs.find(dof_id); if (it == this->dofs.end()) { AKANTU_EXCEPTION("The dof " << dof_id << " does not exists in " << this->id); } return *it->second; } /* -------------------------------------------------------------------------- */ inline const Array & DOFManager::getLocalEquationNumbers(const ID & dof_id) const { return this->getDOFData(dof_id).local_equation_number; } /* -------------------------------------------------------------------------- */ inline Array & DOFManager::getDOFs(const ID & dofs_id) { return *(this->getDOFData(dofs_id).dof); } /* -------------------------------------------------------------------------- */ inline Array & DOFManager::getDOFsDerivatives(const ID & dofs_id, UInt order) { std::vector *> & derivatives = this->getDOFData(dofs_id).dof_derivatives; - if (order >= derivatives.size()) + if (order > derivatives.size()) AKANTU_EXCEPTION("No derivatives of order " << order << " present in " << this->id << " for dof " << dofs_id); - return *derivatives[order]; + return *derivatives[order - 1]; } /* -------------------------------------------------------------------------- */ inline const Array & DOFManager::getSolution(const ID & dofs_id) const { return this->getDOFData(dofs_id).solution; } /* -------------------------------------------------------------------------- */ inline const Array & DOFManager::getBlockedDOFs(const ID & dofs_id) const { return *(this->getDOFData(dofs_id).blocked_dofs); } /* -------------------------------------------------------------------------- */ __END_AKANTU__ #endif /* __AKANTU_DOF_MANAGER_INLINE_IMPL_CC__ */ diff --git a/src/model/integration_scheme/integration_scheme.hh b/src/model/integration_scheme/integration_scheme.hh index 1c69f368a..cdfb5d547 100644 --- a/src/model/integration_scheme/integration_scheme.hh +++ b/src/model/integration_scheme/integration_scheme.hh @@ -1,99 +1,100 @@ /** * @file integration_scheme.hh * * @author Nicolas Richart * * @date Mon Sep 28 10:43:18 2015 * * @brief This class is just a base class for the integration schemes * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_INTEGRATION_SCHEME_HH__ #define __AKANTU_INTEGRATION_SCHEME_HH__ namespace akantu { class DOFManager; } __BEGIN_AKANTU__ class IntegrationScheme { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: enum SolutionType { + _not_defined = -1, _displacement = 0, _temperature = 0, _velocity = 1, _temperature_rate = 1, _acceleration = 2, }; IntegrationScheme(DOFManager & dof_manager, const ID & dof_id, UInt order); virtual ~IntegrationScheme() {} /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// generic interface of a predictor virtual void predictor(Real delta_t) = 0; /// generic interface of a corrector virtual void corrector(const SolutionType & type, Real delta_t) = 0; /// assemble the jacobian matrix virtual void assembleJacobian(const SolutionType & type, Real delta_t) = 0; /// assemble the residual virtual void assembleResidual(bool is_lumped) = 0; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the order of the integration scheme UInt getOrder() const; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// The underlying dofmanager DOFManager & dof_manager; /// The id of the dof treated by this integration scheme. ID dof_id; /// The order of the integrator UInt order; }; __END_AKANTU__ #endif /* __AKANTU_INTEGRATION_SCHEME_HH__ */ diff --git a/src/model/integration_scheme/integration_scheme_1st_order.cc b/src/model/integration_scheme/integration_scheme_1st_order.cc index 139c2b905..99c607d45 100644 --- a/src/model/integration_scheme/integration_scheme_1st_order.cc +++ b/src/model/integration_scheme/integration_scheme_1st_order.cc @@ -1,99 +1,88 @@ /** * @file integration_scheme_1st_order.cc * * @author Nicolas Richart * * @date Fri Oct 23 12:31:32 2015 * * @brief Implementation of the common functions for 1st order time *integrations * * @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 "integration_scheme_1st_order.hh" #include "dof_manager.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ void IntegrationScheme1stOrder::predictor(Real delta_t) { AKANTU_DEBUG_IN(); Array & u = this->dof_manager.getDOFs(this->dof_id); Array & u_dot = this->dof_manager.getDOFsDerivatives(this->dof_id, 1); const Array & blocked_dofs = this->dof_manager.getBlockedDOFs(this->dof_id); this->predictor(delta_t, u, u_dot, blocked_dofs); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void IntegrationScheme1stOrder::corrector(const SolutionType & type, Real delta_t) { AKANTU_DEBUG_IN(); Array & u = this->dof_manager.getDOFs(this->dof_id); Array & u_dot = this->dof_manager.getDOFsDerivatives(this->dof_id, 1); const Array & solution = this->dof_manager.getSolution(this->dof_id); const Array & blocked_dofs = this->dof_manager.getBlockedDOFs(this->dof_id); this->corrector(type, delta_t, u, u_dot, blocked_dofs, solution); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void IntegrationScheme1stOrder::assembleResidual(bool is_lumped) { AKANTU_DEBUG_IN(); const Array & first_derivative = dof_manager.getDOFsDerivatives(this->dof_id, 1); - Array Ctr(first_derivative, true, "Ctr"); - if(!is_lumped) { - const SparseMatrix & C = dof_manager.getMatrix("C"); - Ctr *= C; + if (!is_lumped) { + this->dof_manager.assembleMatMulVectToResidual(this->dof_id, "C", + first_derivative, -1); } else { - const Array & C = dof_manager.getLumpedMatrix("C"); - - UInt nb_dofs = Ctr.getNbComponent() * Ctr.getSize(); - Array::scalar_iterator ctr_it = Ctr.begin_reinterpret(nb_dofs); - Array::scalar_iterator ctr_end = Ctr.end_reinterpret(nb_dofs); - Array::const_scalar_iterator c_it = C.begin_reinterpret(nb_dofs); - - for (; ctr_it != ctr_end; ++ctr_it) { - *ctr_it *= *c_it; - } + this->dof_manager.assembleLumpedMatMulVectToResidual(this->dof_id, "C", + first_derivative, -1); } - dof_manager.assembleToResidual(this->dof_id, Ctr, -1.); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ - __END_AKANTU__ diff --git a/src/model/integration_scheme/integration_scheme_2nd_order.cc b/src/model/integration_scheme/integration_scheme_2nd_order.cc index 044b2bdba..ef01c6df4 100644 --- a/src/model/integration_scheme/integration_scheme_2nd_order.cc +++ b/src/model/integration_scheme/integration_scheme_2nd_order.cc @@ -1,111 +1,101 @@ /** * @file integration_scheme_2nd_order.cc * * @author Nicolas Richart * * @date Tue Oct 20 10:41:12 2015 * * @brief Implementation of the common part of 2nd order integration schemes * * @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 "integration_scheme_2nd_order.hh" #include "dof_manager.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ void IntegrationScheme2ndOrder::predictor(Real delta_t) { AKANTU_DEBUG_IN(); Array & u = this->dof_manager.getDOFs(this->dof_id); Array & u_dot = this->dof_manager.getDOFsDerivatives(this->dof_id, 1); Array & u_dot_dot = this->dof_manager.getDOFsDerivatives(this->dof_id, 2); const Array & blocked_dofs = this->dof_manager.getBlockedDOFs(this->dof_id); this->predictor(delta_t, u, u_dot, u_dot_dot, blocked_dofs); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void IntegrationScheme2ndOrder::corrector(const SolutionType & type, Real delta_t) { AKANTU_DEBUG_IN(); Array & u = this->dof_manager.getDOFs(this->dof_id); Array & u_dot = this->dof_manager.getDOFsDerivatives(this->dof_id, 1); Array & u_dot_dot = this->dof_manager.getDOFsDerivatives(this->dof_id, 2); const Array & solution = this->dof_manager.getSolution(this->dof_id); const Array & blocked_dofs = this->dof_manager.getBlockedDOFs(this->dof_id); this->corrector(type, delta_t, u, u_dot, u_dot_dot, blocked_dofs, solution); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void IntegrationScheme2ndOrder::assembleResidual(bool is_lumped) { AKANTU_DEBUG_IN(); try { const Array & first_derivative = - dof_manager.getDOFsDerivatives(this->dof_id, 1); - Array Cv(first_derivative, true, "Cv"); - const SparseMatrix & C = dof_manager.getMatrix("C"); - Cv *= C; - dof_manager.assembleToResidual(this->dof_id, Cv, -1.); + this->dof_manager.getDOFsDerivatives(this->dof_id, 1); + + this->dof_manager.assembleMatMulVectToResidual(this->dof_id, "C", + first_derivative, -1); } catch (...) { } const Array & second_derivative = - dof_manager.getDOFsDerivatives(this->dof_id, 2); - Array Ma(second_derivative, true, "Ma"); + this->dof_manager.getDOFsDerivatives(this->dof_id, 2); + if (!is_lumped) { - const SparseMatrix & M = dof_manager.getMatrix("M"); - Ma *= M; + this->dof_manager.assembleMatMulVectToResidual(this->dof_id, "M", + second_derivative, -1); } else { - const Array & M = dof_manager.getLumpedMatrix("M"); - - UInt nb_dofs = Ma.getNbComponent() * Ma.getSize(); - Array::scalar_iterator ma_it = Ma.begin_reinterpret(nb_dofs); - Array::scalar_iterator ma_end = Ma.end_reinterpret(nb_dofs); - Array::const_scalar_iterator m_it = M.begin_reinterpret(nb_dofs); - - for (; ma_it != ma_end; ++ma_it) { - *ma_it *= *m_it; - } + this->dof_manager.assembleLumpedMatMulVectToResidual(this->dof_id, "M", + second_derivative, -1); } - dof_manager.assembleToResidual(this->dof_id, Ma, -1.); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/model_solver.cc b/src/model/model_solver.cc index c1ffa4724..a2d0866a6 100644 --- a/src/model/model_solver.cc +++ b/src/model/model_solver.cc @@ -1,470 +1,503 @@ /** * @file model_solver.cc * * @author Nicolas Richart * * @date Wed Aug 12 13:31:56 2015 * * @brief Implementation of ModelSolver * * @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 "model_solver.hh" #include "mesh.hh" #include "dof_manager.hh" #if defined(AKANTU_USE_MPI) #include "mpi_type_wrapper.hh" #endif #if defined(AKANTU_USE_MUMPS) #include "dof_manager_default.hh" #endif #if defined(AKANTU_USE_PETSC) #include "dof_manager_petsc.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ ModelSolver::ModelSolver(Mesh & mesh, const ID & id, UInt memory_id) : Parsable(_st_solver, id), SolverCallback(), parent_id(id), - parent_memory_id(memory_id), mesh(mesh), dof_manager(NULL) {} + parent_memory_id(memory_id), mesh(mesh), dof_manager(NULL), + default_solver_id("") {} /* -------------------------------------------------------------------------- */ ModelSolver::~ModelSolver() { delete this->dof_manager; } /* -------------------------------------------------------------------------- */ void ModelSolver::initDOFManager() { std::pair sub_sect = getStaticParser().getSubSections(_st_solver); if (sub_sect.first != sub_sect.second) { AKANTU_EXCEPTION("More than on solver section present in the input file"); } const ParserSection & section = *sub_sect.first; std::string solver_type = section.getName(); this->initDOFManager(solver_type); } /* -------------------------------------------------------------------------- */ void ModelSolver::initDOFManager(const ID & solver_type) { if (solver_type == "petsc") { #if defined(AKANTU_USE_PETSC) ID id = this->parent_id + ":dof_manager_petsc"; this->dof_manager = new DOFManagerPETSc(id, this->parent_memory_id); #else AKANTU_EXCEPTION( "To use PETSc you have to activate it in the compilations options"); #endif } else if (solver_type == "mumps") { #if defined(AKANTU_USE_MUMPS) ID id = this->parent_id + ":dof_manager_default"; this->dof_manager = new DOFManagerDefault(id, this->parent_memory_id); #else AKANTU_EXCEPTION( "To use MUMPS you have to activate it in the compilations options"); #endif } else { AKANTU_EXCEPTION( "To use the solver " << solver_type << " you will have to code it. This is an unknown solver type."); } this->dof_manager->registerMesh(mesh); this->setDOFManager(*this->dof_manager); } /* -------------------------------------------------------------------------- */ -void ModelSolver::solveStep(ID solver_id) { - AKANTU_DEBUG_IN(); +TimeStepSolver & ModelSolver::getSolver(const ID & solver_id) { + ID tmp_solver_id = solver_id; + if (tmp_solver_id == "") + tmp_solver_id = this->default_solver_id; + + TimeStepSolver & tss = this->dof_manager->getTimeStepSolver(tmp_solver_id); + return tss; +} +/* -------------------------------------------------------------------------- */ +const TimeStepSolver & ModelSolver::getSolver(const ID & solver_id) const { + ID tmp_solver_id = solver_id; if (solver_id == "") - solver_id = default_solver_id; + tmp_solver_id = this->default_solver_id; - TimeStepSolver & tss = this->dof_manager->getTimeStepSolver(solver_id); + const TimeStepSolver & tss = + this->dof_manager->getTimeStepSolver(tmp_solver_id); + return tss; +} + +/* -------------------------------------------------------------------------- */ +void ModelSolver::solveStep(const ID & solver_id) { + AKANTU_DEBUG_IN(); + TimeStepSolver & tss = this->getSolver(solver_id); // make one non linear solve tss.solveStep(*this); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ModelSolver::getNewSolver(const ID & solver_id, TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) { if (this->default_solver_id == "") { this->default_solver_id = solver_id; } if (non_linear_solver_type == _nls_auto) { switch (time_step_solver_type) { case _tsst_dynamic: case _tsst_static: non_linear_solver_type = _nls_newton_raphson; break; case _tsst_dynamic_lumped: non_linear_solver_type = _nls_lumped; break; } } NonLinearSolver & nls = this->dof_manager->getNewNonLinearSolver( solver_id, non_linear_solver_type); this->dof_manager->getNewTimeStepSolver(solver_id, time_step_solver_type, nls); } +/* -------------------------------------------------------------------------- */ +Real ModelSolver::getTimeStep(const ID & solver_id) const { + const TimeStepSolver & tss = this->getSolver(solver_id); + + return tss.getTimeStep(); +} + +/* -------------------------------------------------------------------------- */ +void ModelSolver::setTimeStep(Real time_step, const ID & solver_id) { + TimeStepSolver & tss = this->getSolver(solver_id); + + return tss.setTimeStep(time_step); +} + /* -------------------------------------------------------------------------- */ void ModelSolver::setIntegrationScheme( const ID & solver_id, const ID & dof_id, - const IntegrationSchemeType & integration_scheme_type) { + const IntegrationSchemeType & integration_scheme_type, + IntegrationScheme::SolutionType solution_type) { TimeStepSolver & tss = this->dof_manager->getTimeStepSolver(solver_id); - tss.setIntegrationScheme(dof_id, integration_scheme_type); + tss.setIntegrationScheme(dof_id, integration_scheme_type, solution_type); } /* -------------------------------------------------------------------------- */ void ModelSolver::predictor() {} /* -------------------------------------------------------------------------- */ void ModelSolver::corrector() {} /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ // void ModelSolver::setIntegrationScheme( // const ID & solver_id, const ID & dof_id, // const IntegrationSchemeType & integration_scheme_type) {} /* -------------------------------------------------------------------------- */ // /* -------------------------------------------------------------------------- // */ // template // void SolidMechanicsModel::solve(Array &increment, Real block_val, // bool need_factorize, bool // has_profile_changed) { // if(has_profile_changed) { // this->initJacobianMatrix(); // need_factorize = true; // } // updateResidualInternal(); //doesn't do anything for static // if(need_factorize) { // Real c = 0.,d = 0.,e = 0.; // if(method == _static) { // AKANTU_DEBUG_INFO("Solving K inc = r"); // e = 1.; // } else { // AKANTU_DEBUG_INFO("Solving (c M + d C + e K) inc = r"); // NewmarkBeta * nmb_int = dynamic_cast(integrator); // c = nmb_int->getAccelerationCoefficient(time_step); // d = nmb_int->getVelocityCoefficient(time_step); // e = nmb_int->getDisplacementCoefficient(time_step); // } // jacobian_matrix->clear(); // // J = c M + d C + e K // if(stiffness_matrix) // jacobian_matrix->add(*stiffness_matrix, e); // if(mass_matrix) // jacobian_matrix->add(*mass_matrix, c); // #if !defined(AKANTU_NDEBUG) // if(mass_matrix && AKANTU_DEBUG_TEST(dblDump)) // mass_matrix->saveMatrix("M.mtx"); // #endif // if(velocity_damping_matrix) // jacobian_matrix->add(*velocity_damping_matrix, d); // jacobian_matrix->applyBoundary(*blocked_dofs, block_val); // #if !defined(AKANTU_NDEBUG) // if(AKANTU_DEBUG_TEST(dblDump)) // jacobian_matrix->saveMatrix("J.mtx"); // #endif // solver->factorize(); // } // // if (rhs.getSize() != 0) // // solver->setRHS(rhs); // // else // solver->setOperators(); // solver->setRHS(*residual); // // solve @f[ J \delta w = r @f] // solver->solve(increment); // UInt nb_nodes = displacement-> getSize(); // UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes; // bool * blocked_dofs_val = blocked_dofs->storage(); // Real * increment_val = increment.storage(); // for (UInt j = 0; j < nb_degree_of_freedom; // ++j,++increment_val, ++blocked_dofs_val) { // if ((*blocked_dofs_val)) // *increment_val = 0.0; // } // } // /* -------------------------------------------------------------------------- // */ // template // bool SolidMechanicsModel::solveStatic(Real tolerance, UInt max_iteration, // bool do_not_factorize) { // AKANTU_DEBUG_INFO("Solving Ku = f"); // AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL, // "You should first initialize the implicit solver and // assemble the stiffness matrix by calling // initImplicit"); // AnalysisMethod analysis_method=method; // Real error = 0.; // method=_static; // bool converged = this->template solveStep(tolerance, // error, max_iteration, do_not_factorize); // method=analysis_method; // return converged; // } // /* -------------------------------------------------------------------------- // */ // template // bool SolidMechanicsModel::solveStep(Real tolerance, // UInt max_iteration) { // Real error = 0.; // return this->template solveStep(tolerance, // error, // max_iteration); // } // /* -------------------------------------------------------------------------- // */ // template // bool SolidMechanicsModel::solveStep(Real tolerance, Real & error, UInt // max_iteration, // bool do_not_factorize) { // EventManager::sendEvent(SolidMechanicsModelEvent::BeforeSolveStepEvent(method)); // this->implicitPred(); // this->updateResidual(); // AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL, // "You should first initialize the implicit solver and // assemble the stiffness matrix"); // bool need_factorize = !do_not_factorize; // if (method==_implicit_dynamic) { // AKANTU_DEBUG_ASSERT(mass_matrix != NULL, // "You should first initialize the implicit solver and // assemble the mass matrix"); // } // switch (cmethod) { // case _scm_newton_raphson_tangent: // case _scm_newton_raphson_tangent_not_computed: // break; // case _scm_newton_raphson_tangent_modified: // this->assembleStiffnessMatrix(); // break; // default: // AKANTU_DEBUG_ERROR("The resolution method " << cmethod << " has not been // implemented!"); // } // this->n_iter = 0; // bool converged = false; // error = 0.; // if(criteria == _scc_residual) { // converged = this->testConvergence (tolerance, error); // if(converged) return converged; // } // do { // if (cmethod == _scm_newton_raphson_tangent) // this->assembleStiffnessMatrix(); // solve (*increment, 1., // need_factorize); // this->implicitCorr(); // if(criteria == _scc_residual) this->updateResidual(); // converged = this->testConvergence (tolerance, error); // if(criteria == _scc_increment && !converged) this->updateResidual(); // //this->dump(); // this->n_iter++; // AKANTU_DEBUG_INFO("[" << criteria << "] Convergence iteration " // << std::setw(std::log10(max_iteration)) << this->n_iter // << ": error " << error << (converged ? " < " : " > ") // << tolerance); // switch (cmethod) { // case _scm_newton_raphson_tangent: // need_factorize = true; // break; // case _scm_newton_raphson_tangent_not_computed: // case _scm_newton_raphson_tangent_modified: // need_factorize = false; // break; // default: // AKANTU_DEBUG_ERROR("The resolution method " << cmethod << " has not // been implemented!"); // } // } while (!converged && this->n_iter < max_iteration); // // this makes sure that you have correct strains and stresses after the // solveStep function (e.g., for dumping) // if(criteria == _scc_increment) this->updateResidual(); // if (converged) { // EventManager::sendEvent(SolidMechanicsModelEvent::AfterSolveStepEvent(method)); // } else if(this->n_iter == max_iteration) { // AKANTU_DEBUG_WARNING("[" << criteria << "] Convergence not reached after // " // << std::setw(std::log10(max_iteration)) << // this->n_iter << // " iteration" << (this->n_iter == 1 ? "" : "s") << // "!" << std::endl); // } // return converged; // } // void SolidMechanicsModel::updateResidualInternal() { // AKANTU_DEBUG_IN(); // AKANTU_DEBUG_INFO("Update the residual"); // // f = f_ext - f_int - Ma - Cv = r - Ma - Cv; // if(method != _static) { // // f -= Ma // if(mass_matrix) { // // if full mass_matrix // Array * Ma = new Array(*acceleration, true, "Ma"); // *Ma *= *mass_matrix; // /// \todo check unit conversion for implicit dynamics // // *Ma /= f_m2a // *residual -= *Ma; // delete Ma; // } else if (mass) { // // else lumped mass // UInt nb_nodes = acceleration->getSize(); // UInt nb_degree_of_freedom = acceleration->getNbComponent(); // Real * mass_val = mass->storage(); // Real * accel_val = acceleration->storage(); // Real * res_val = residual->storage(); // bool * blocked_dofs_val = blocked_dofs->storage(); // for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) { // if(!(*blocked_dofs_val)) { // *res_val -= *accel_val * *mass_val /f_m2a; // } // blocked_dofs_val++; // res_val++; // mass_val++; // accel_val++; // } // } else { // AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix."); // } // // f -= Cv // if(velocity_damping_matrix) { // Array * Cv = new Array(*velocity); // *Cv *= *velocity_damping_matrix; // *residual -= *Cv; // delete Cv; // } // } // AKANTU_DEBUG_OUT(); // } // /* -------------------------------------------------------------------------- // */ // void SolidMechanicsModel::solveLumped(Array & x, // const Array & A, // const Array & b, // const Array & blocked_dofs, // Real alpha) { // Real * A_val = A.storage(); // Real * b_val = b.storage(); // Real * x_val = x.storage(); // bool * blocked_dofs_val = blocked_dofs.storage(); // UInt nb_degrees_of_freedom = x.getSize() * x.getNbComponent(); // for (UInt n = 0; n < nb_degrees_of_freedom; ++n) { // if(!(*blocked_dofs_val)) { // *x_val = alpha * (*b_val / *A_val); // } // x_val++; // A_val++; // b_val++; // blocked_dofs_val++; // } // } /* -------------------------------------------------------------------------- */ // void TimeStepSolverDefault::updateAcceleration() { // AKANTU_DEBUG_IN(); // updateResidualInternal(); // if (method == _explicit_lumped_mass) { // /* residual = residual_{n+1} - M * acceleration_n therefore // solution = increment acceleration not acceleration */ // solveLumped(*increment_acceleration, *mass, *residual, *blocked_dofs, // f_m2a); // } else if (method == _explicit_consistent_mass) { // solve(*increment_acceleration); // } // AKANTU_DEBUG_OUT(); // } __END_AKANTU__ diff --git a/src/model/model_solver.hh b/src/model/model_solver.hh index b63c0b787..6a08ca926 100644 --- a/src/model/model_solver.hh +++ b/src/model/model_solver.hh @@ -1,126 +1,139 @@ /** * @file model_solver.hh * * @author Nicolas Richart * * @date Wed Jul 22 10:53:10 2015 * * @brief Class regrouping the common solve interface to the different models * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "parsable.hh" #include "solver_callback.hh" +#include "integration_scheme.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MODEL_SOLVER_HH__ #define __AKANTU_MODEL_SOLVER_HH__ namespace akantu { - class Mesh; - class DOFManager; - class IntegrationScheme; +class Mesh; +class DOFManager; +class TimeStepSolver; } __BEGIN_AKANTU__ class ModelSolver : public Parsable, public SolverCallback { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: ModelSolver(Mesh & mesh, const ID & id, UInt memory_id); virtual ~ModelSolver(); /// initialize the dof manager based on solver type passed in the input file void initDOFManager(); /// initialize the dof manager based on the used chosen solver type void initDOFManager(const ID & solver_type); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// solve a step using a given pre instantiated time step solver and /// nondynamic linear solver - void solveStep(ID time_step_solver_id = ""); + void solveStep(const ID & solver_id = ""); /// Initialize a time solver that can be used afterwards with its id - void - getNewSolver(const ID & solver_id, - TimeStepSolverType time_step_solver_type, - NonLinearSolverType non_linear_solver_type = _nls_auto); + void getNewSolver(const ID & solver_id, + TimeStepSolverType time_step_solver_type, + NonLinearSolverType non_linear_solver_type = _nls_auto); /// set an integration scheme for a given dof and a given solver void setIntegrationScheme(const ID & solver_id, const ID & dof_id, - const IntegrationSchemeType & integration_scheme_type); + const IntegrationSchemeType & integration_scheme_type, + IntegrationScheme::SolutionType solution_type = + IntegrationScheme::_not_defined); /// set an externally instantiated integration scheme void setIntegrationScheme(const ID & solver_id, const ID & dof_id, - IntegrationScheme & integration_scheme); + IntegrationScheme & integration_scheme, + IntegrationScheme::SolutionType solution_type = + IntegrationScheme::_not_defined); /* ------------------------------------------------------------------------ */ /* SolverCallback interface */ /* ------------------------------------------------------------------------ */ public: /// Predictor interface for the callback virtual void predictor(); /// Corrector interface for the callback virtual void corrector(); /// AssembleResidual interface for the callback virtual void assembleResidual() = 0; /// AssembleJacobian interface for the callback virtual void assembleJacobian() = 0; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: DOFManager & getDOFManager() { return *this->dof_manager; } + /// get the time step of a given solver + Real getTimeStep(const ID & solver_id = "") const; + /// set the time step of a given solver + void setTimeStep(Real time_step, const ID & solver_id = ""); + +private: + TimeStepSolver & getSolver(const ID & solver_id); + const TimeStepSolver & getSolver(const ID & solver_id) const; + /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: ID parent_id; UInt parent_memory_id; /// Underlying mesh Mesh & mesh; /// Underlying dof_manager (the brain...) DOFManager * dof_manager; /// Default time step solver to use ID default_solver_id; }; __END_AKANTU__ #endif /* __AKANTU_MODEL_SOLVER_HH__ */ diff --git a/src/model/non_linear_solver_newton_raphson.cc b/src/model/non_linear_solver_newton_raphson.cc index 0a2f691fb..697bdb6f0 100644 --- a/src/model/non_linear_solver_newton_raphson.cc +++ b/src/model/non_linear_solver_newton_raphson.cc @@ -1,161 +1,163 @@ /** * @file non_linear_solver_newton_raphson.cc * * @author Nicolas Richart * * @date Tue Aug 25 00:57:00 2015 * * @brief Implementation of the default NonLinearSolver * * @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 "non_linear_solver_newton_raphson.hh" #include "dof_manager_default.hh" #include "solver_callback.hh" #include "static_communicator.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ NonLinearSolverNewtonRaphson::NonLinearSolverNewtonRaphson( DOFManagerDefault & dof_manager, const NonLinearSolverType & non_linear_solver_type, const ID & id, UInt memory_id) : NonLinearSolver(dof_manager, non_linear_solver_type, id, memory_id), dof_manager(dof_manager), - solver(dof_manager, "jacobian", id + ":sparse_solver", memory_id), + solver(dof_manager, "J", id + ":sparse_solver", memory_id), convergence_criteria_type(_scc_solution), convergence_criteria(1e-10), max_iterations(10), n_iter(0), error(0.), converged(false) { this->supported_type.insert(_nls_newton_raphson_modified); this->supported_type.insert(_nls_newton_raphson); this->supported_type.insert(_nls_linear); this->checkIfTypeIsSupported(); } /* -------------------------------------------------------------------------- */ NonLinearSolverNewtonRaphson::~NonLinearSolverNewtonRaphson() {} /* ------------------------------------------------------------------------ */ void NonLinearSolverNewtonRaphson::solve(SolverCallback & solver_callback) { solver_callback.predictor(); + solver_callback.assembleResidual(); + if(this->non_linear_solver_type == _nls_newton_raphson_modified) solver_callback.assembleJacobian(); this->n_iter = 0; this->converged = false; + if (this->convergence_criteria_type == _scc_residual) { this->converged = this->testConvergence(this->dof_manager.getResidual()); if (this->converged) return; } do { if (this->non_linear_solver_type == _nls_newton_raphson) solver_callback.assembleJacobian(); this->solver.solve(); solver_callback.corrector(); // EventManager::sendEvent(NonLinearSolver::AfterSparseSolve(method)); if (this->convergence_criteria_type == _scc_residual) { solver_callback.assembleResidual(); this->converged = this->testConvergence(this->dof_manager.getResidual()); } else { this->converged = this->testConvergence(this->dof_manager.getGlobalSolution()); } if (this->convergence_criteria_type == _scc_solution && !this->converged) solver_callback.assembleResidual(); // this->dump(); this->n_iter++; AKANTU_DEBUG_INFO( "[" << this->convergence_criteria_type << "] Convergence iteration " << std::setw(std::log10(this->max_iterations)) << this->n_iter << ": error " << this->error << (this->converged ? " < " : " > ") << this->convergence_criteria); } while (!this->converged && this->n_iter < this->max_iterations); // this makes sure that you have correct strains and stresses after the // solveStep function (e.g., for dumping) if (this->convergence_criteria_type == _scc_solution) solver_callback.assembleResidual(); if (this->converged) { // EventManager::sendEvent( // SolidMechanicsModelEvent::AfterNonLinearSolverSolves(method)); } else if (this->n_iter == this->max_iterations) { AKANTU_DEBUG_WARNING("[" << this->convergence_criteria_type << "] Convergence not reached after " << std::setw(std::log10(this->max_iterations)) << this->n_iter << " iteration" - << (this->n_iter == 1 ? "" : "s") << "!" - << std::endl); + << (this->n_iter == 1 ? "" : "s") << "!"); } return; } /* -------------------------------------------------------------------------- */ bool NonLinearSolverNewtonRaphson::testConvergence(const Array & array) { AKANTU_DEBUG_IN(); const Array & blocked_dofs = this->dof_manager.getGlobalBlockedDOFs(); UInt nb_degree_of_freedoms = array.getSize(); Array::const_scalar_iterator arr_it = array.begin(); Array::const_scalar_iterator bld_it = blocked_dofs.begin(); Real norm = 0.; for (UInt n = 0; n < nb_degree_of_freedoms; ++n, ++arr_it, ++bld_it) { bool is_local_node = this->dof_manager.isLocalOrMasterDOF(n); if (!(*bld_it) && is_local_node) { norm += *arr_it * *arr_it; } } StaticCommunicator::getStaticCommunicator().allReduce(norm, _so_sum); norm = std::sqrt(norm); AKANTU_DEBUG_ASSERT(!Math::isnan(norm), - "Something goes wrong in the solve phase"); + "Something went wrong in the solve phase"); this->error = norm; return (error < this->convergence_criteria); } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/time_step_solver.hh b/src/model/time_step_solver.hh index aed6aded9..c0a5ecd1a 100644 --- a/src/model/time_step_solver.hh +++ b/src/model/time_step_solver.hh @@ -1,103 +1,112 @@ /** * @file time_step_solver.hh * * @author Nicolas Richart * * @date Mon Aug 24 12:42:04 2015 * * @brief This corresponding to the time step evolution solver * * @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_memory.hh" #include "aka_array.hh" #include "solver_callback.hh" - +#include "integration_scheme.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_TIME_STEP_SOLVER_HH__ #define __AKANTU_TIME_STEP_SOLVER_HH__ namespace akantu { class DOFManager; class NonLinearSolver; } __BEGIN_AKANTU__ class TimeStepSolver : public Memory, public SolverCallback { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: TimeStepSolver(DOFManager & dof_manager, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver, const ID & id, UInt memory_id); virtual ~TimeStepSolver(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// solves on step virtual void solveStep(SolverCallback & solver_callback) = 0; /// register an integration scheme for a given dof - virtual void setIntegrationScheme(const ID & dof_id, - const IntegrationSchemeType & type) = 0; + virtual void + setIntegrationScheme(const ID & dof_id, const IntegrationSchemeType & type, + IntegrationScheme::SolutionType + solution_type = IntegrationScheme::_not_defined) = 0; /* ------------------------------------------------------------------------ */ /* Solver Callback interface */ /* ------------------------------------------------------------------------ */ public: /// implementation of the SolverCallback::predictor() virtual void predictor(); /// implementation of the SolverCallback::corrector() virtual void corrector(); /// implementation of the SolverCallback::assembleJacobian() virtual void assembleJacobian(); /// implementation of the SolverCallback::assembleResidual() virtual void assembleResidual(); + /* ------------------------------------------------------------------------ */ + /* Accessor */ + /* ------------------------------------------------------------------------ */ +public: + AKANTU_GET_MACRO(TimeStep, time_step, Real); + AKANTU_SET_MACRO(TimeStep, time_step, Real); + /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// Underlying dof manager containing the dof to treat DOFManager & _dof_manager; /// Type of solver TimeStepSolverType type; /// The time step for this solver Real time_step; /// Temporary storage for solver callback SolverCallback * solver_callback; /// NonLinearSolver used by this tome step solver NonLinearSolver & non_linear_solver; }; __END_AKANTU__ #endif /* __AKANTU_TIME_STEP_SOLVER_HH__ */ diff --git a/src/model/time_step_solvers/time_step_solver.cc b/src/model/time_step_solvers/time_step_solver.cc index fcfdbb7cc..e23a7f0eb 100644 --- a/src/model/time_step_solvers/time_step_solver.cc +++ b/src/model/time_step_solvers/time_step_solver.cc @@ -1,88 +1,91 @@ /** * @file time_step_solver.cc * * @author Nicolas Richart * * @date Mon Oct 12 16:56:43 2015 * * @brief Implementation of common part of TimeStepSolvers * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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 "time_step_solver.hh" #include "non_linear_solver.hh" +#include "dof_manager.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ TimeStepSolver::TimeStepSolver(DOFManager & dof_manager, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver, const ID & id, UInt memory_id) : Memory(id, memory_id), SolverCallback(dof_manager), - _dof_manager(dof_manager), type(type), solver_callback(NULL), + _dof_manager(dof_manager), type(type), time_step(0.), solver_callback(NULL), non_linear_solver(non_linear_solver) {} /* -------------------------------------------------------------------------- */ TimeStepSolver::~TimeStepSolver() {} /* -------------------------------------------------------------------------- */ void TimeStepSolver::predictor() { AKANTU_DEBUG_ASSERT( this->solver_callback != NULL, "This function cannot be called if the solver_callback is not set"); this->solver_callback->predictor(); } /* -------------------------------------------------------------------------- */ void TimeStepSolver::corrector() { AKANTU_DEBUG_ASSERT( this->solver_callback != NULL, "This function cannot be called if the solver_callback is not set"); this->solver_callback->corrector(); } /* -------------------------------------------------------------------------- */ void TimeStepSolver::assembleJacobian() { AKANTU_DEBUG_ASSERT( this->solver_callback != NULL, "This function cannot be called if the solver_callback is not set"); + this->_dof_manager.clearJacobian(); this->solver_callback->assembleJacobian(); } /* -------------------------------------------------------------------------- */ void TimeStepSolver::assembleResidual() { AKANTU_DEBUG_ASSERT( this->solver_callback != NULL, "This function cannot be called if the solver_callback is not set"); + this->_dof_manager.clearResidual(); this->solver_callback->assembleResidual(); } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/time_step_solvers/time_step_solver_default.cc b/src/model/time_step_solvers/time_step_solver_default.cc index 138f98f8f..8ccacec37 100644 --- a/src/model/time_step_solvers/time_step_solver_default.cc +++ b/src/model/time_step_solvers/time_step_solver_default.cc @@ -1,248 +1,257 @@ /** * @file time_step_solver_default.cc * * @author Nicolas Richart * * @date Wed Sep 16 10:20:55 2015 * * @brief Default implementation of the time step solver * * @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 "time_step_solver_default.hh" #include "mesh.hh" #include "dof_manager_default.hh" #include "sparse_matrix_aij.hh" #include "pseudo_time.hh" #include "integration_scheme_1st_order.hh" #include "integration_scheme_2nd_order.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ TimeStepSolverDefault::TimeStepSolverDefault( DOFManagerDefault & dof_manager, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver, const ID & id, UInt memory_id) : TimeStepSolver(dof_manager, type, non_linear_solver, id, memory_id), dof_manager(dof_manager), is_mass_lumped(false) { switch (type) { case _tsst_dynamic: break; case _tsst_dynamic_lumped: this->is_mass_lumped = true; break; case _tsst_static: /// initialize a static time solver for allback dofs break; } } /* -------------------------------------------------------------------------- */ void TimeStepSolverDefault::setIntegrationScheme( - const ID & dof_id, const IntegrationSchemeType & type) { + const ID & dof_id, const IntegrationSchemeType & type, + IntegrationScheme::SolutionType solution_type) { if (this->integration_schemes.find(dof_id) != this->integration_schemes.end()) { AKANTU_EXCEPTION("Their DOFs " << dof_id << " have already an integration scheme associated"); } IntegrationScheme * integration_scheme = NULL; if (this->is_mass_lumped) { switch (type) { case _ist_forward_euler: { integration_scheme = new ForwardEuler(dof_manager, dof_id); break; } case _ist_central_difference: { integration_scheme = new CentralDifference(dof_manager, dof_id); break; } default: AKANTU_EXCEPTION( "This integration scheme cannot be used in lumped dynamic"); } } else { switch (type) { case _ist_pseudo_time: { integration_scheme = new PseudoTime(dof_manager, dof_id); break; } case _ist_forward_euler: { integration_scheme = new ForwardEuler(dof_manager, dof_id); break; } case _ist_trapezoidal_rule_1: { integration_scheme = new TrapezoidalRule1(dof_manager, dof_id); break; } case _ist_backward_euler: { integration_scheme = new BackwardEuler(dof_manager, dof_id); break; } case _ist_central_difference: { integration_scheme = new CentralDifference(dof_manager, dof_id); break; } case _ist_fox_goodwin: { integration_scheme = new FoxGoodwin(dof_manager, dof_id); break; } case _ist_trapezoidal_rule_2: { integration_scheme = new TrapezoidalRule2(dof_manager, dof_id); break; } case _ist_linear_acceleration: { integration_scheme = new LinearAceleration(dof_manager, dof_id); break; } // Write a c++11 version of the constructor with initializer list that // contains the arguments for the integration scheme case _ist_generalized_trapezoidal: case _ist_newmark_beta: AKANTU_EXCEPTION( "This time step solvers cannot be created with this constructor"); } } this->integration_schemes[dof_id] = integration_scheme; + this->solution_types[dof_id] = solution_type; + this->integration_schemes_owner.insert(dof_id); } /* -------------------------------------------------------------------------- */ TimeStepSolverDefault::~TimeStepSolverDefault() { DOFsIntegrationSchemesOwner::iterator it = this->integration_schemes_owner.begin(); DOFsIntegrationSchemesOwner::iterator end = this->integration_schemes_owner.end(); for (; it != end; ++it) { delete this->integration_schemes[*it]; } this->integration_schemes.clear(); } /* -------------------------------------------------------------------------- */ void TimeStepSolverDefault::solveStep(SolverCallback & solver_callback) { this->solver_callback = &solver_callback; this->non_linear_solver.solve(*this); this->solver_callback = NULL; } /* -------------------------------------------------------------------------- */ void TimeStepSolverDefault::predictor() { AKANTU_DEBUG_IN(); TimeStepSolver::predictor(); DOFsIntegrationSchemes::iterator integration_scheme_it = this->integration_schemes.begin(); DOFsIntegrationSchemes::iterator integration_scheme_end = this->integration_schemes.end(); for (; integration_scheme_it != integration_scheme_end; ++integration_scheme_it) { integration_scheme_it->second->predictor(this->time_step); } // UInt nb_degree_of_freedom = u.getSize() * u.getNbComponent(); // Array::scalar_iterator incr_it = // increment.begin_reinterpret(nb_degree_of_freedom); // Array::const_scalar_iterator u_it = // u.begin_reinterpret(nb_degree_of_freedom); // Array::const_scalar_iterator u_end = // u.end_reinterpret(nb_degree_of_freedom); // for (; u_it != u_end; ++u_it, ++incr_it) { // *incr_it = *u_it - *incr_it; // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void TimeStepSolverDefault::corrector() { AKANTU_DEBUG_IN(); TimeStepSolver::corrector(); DOFsIntegrationSchemes::iterator integration_scheme_it = this->integration_schemes.begin(); DOFsIntegrationSchemes::iterator integration_scheme_end = this->integration_schemes.end(); for (; integration_scheme_it != integration_scheme_end; ++integration_scheme_it) { + IntegrationScheme::SolutionType solution_type = + this->solution_types[integration_scheme_it->first]; + integration_scheme_it->second->corrector( - IntegrationScheme::SolutionType(this->solution_type), this->time_step); + solution_type, this->time_step); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void TimeStepSolverDefault::assembleJacobian() { AKANTU_DEBUG_IN(); TimeStepSolver::assembleJacobian(); DOFsIntegrationSchemes::iterator integration_scheme_it = this->integration_schemes.begin(); DOFsIntegrationSchemes::iterator integration_scheme_end = this->integration_schemes.end(); for (; integration_scheme_it != integration_scheme_end; ++integration_scheme_it) { - integration_scheme_it->second->assembleJacobian( - IntegrationScheme::SolutionType(this->solution_type), this->time_step); + IntegrationScheme::SolutionType solution_type = + this->solution_types[integration_scheme_it->first]; + + integration_scheme_it->second->assembleJacobian(solution_type, + this->time_step); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void TimeStepSolverDefault::assembleResidual() { AKANTU_DEBUG_IN(); TimeStepSolver::assembleResidual(); DOFsIntegrationSchemes::iterator integration_scheme_it = this->integration_schemes.begin(); DOFsIntegrationSchemes::iterator integration_scheme_end = this->integration_schemes.end(); for (; integration_scheme_it != integration_scheme_end; ++integration_scheme_it) { integration_scheme_it->second->assembleResidual(this->is_mass_lumped); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/time_step_solvers/time_step_solver_default.hh b/src/model/time_step_solvers/time_step_solver_default.hh index feab3a9d8..6acccca44 100644 --- a/src/model/time_step_solvers/time_step_solver_default.hh +++ b/src/model/time_step_solvers/time_step_solver_default.hh @@ -1,105 +1,109 @@ /** * @file time_step_solver_default.hh * * @author Nicolas Richart * * @date Mon Aug 24 17:10:29 2015 * * @brief Default implementation for the time stepper * * @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 "time_step_solver.hh" +#include "integration_scheme.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_TIME_STEP_SOLVER_DEFAULT_HH__ #define __AKANTU_TIME_STEP_SOLVER_DEFAULT_HH__ namespace akantu { -class IntegrationScheme; class DOFManagerDefault; } __BEGIN_AKANTU__ class TimeStepSolverDefault : public TimeStepSolver { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: TimeStepSolverDefault(DOFManagerDefault & dof_manager, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver, const ID & id, UInt memory_id); virtual ~TimeStepSolverDefault(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// registers an integration scheme for a given dof void setIntegrationScheme(const ID & dof_id, - const IntegrationSchemeType & type); + const IntegrationSchemeType & type, + IntegrationScheme::SolutionType solution_type = + IntegrationScheme::_not_defined); /// implementation of the TimeStepSolver::predictor() virtual void predictor(); /// implementation of the TimeStepSolver::corrector() virtual void corrector(); /// implementation of the TimeStepSolver::assembleJacobian() virtual void assembleJacobian(); /// implementation of the TimeStepSolver::assembleResidual() virtual void assembleResidual(); /// implementation of the generic TimeStepSolver::solveStep() virtual void solveStep(SolverCallback & solver_callback); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: typedef std::map DOFsIntegrationSchemes; + typedef std::map + DOFsIntegrationSchemesSolutionTypes; typedef std::set DOFsIntegrationSchemesOwner; /// DOFManager with its real type DOFManagerDefault & dof_manager; /// Underlying integration scheme per dof, \todo check what happens in dynamic /// in case of coupled equations DOFsIntegrationSchemes integration_schemes; /// defines if the solver is owner of the memory or not DOFsIntegrationSchemesOwner integration_schemes_owner; /// Type of corrector to use - UInt solution_type; + DOFsIntegrationSchemesSolutionTypes solution_types; /// define if the mass matrix is lumped or not bool is_mass_lumped; }; __END_AKANTU__ #endif /* __AKANTU_TIME_STEP_SOLVER_DEFAULT_HH__ */ diff --git a/src/solver/solver_mumps.cc b/src/solver/solver_mumps.cc index 4a6f7a091..e9b237652 100644 --- a/src/solver/solver_mumps.cc +++ b/src/solver/solver_mumps.cc @@ -1,403 +1,407 @@ /** * @file sparse_solver_mumps.cc * * @author Nicolas Richart * * @date creation: Mon Dec 13 2010 * @date last modification: Tue Jan 19 2016 * * @brief implem of SparseSolverMumps class * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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 * * @subsection Ctrl_param Control parameters * * ICNTL(1), * ICNTL(2), * ICNTL(3) : output streams for error, diagnostics, and global messages * * ICNTL(4) : verbose level : 0 no message - 4 all messages * * ICNTL(5) : type of matrix, 0 assembled, 1 elementary * * ICNTL(6) : control the permutation and scaling(default 7) see mumps doc for * more information * * ICNTL(7) : determine the pivot order (default 7) see mumps doc for more * information * * ICNTL(8) : describe the scaling method used * * ICNTL(9) : 1 solve A x = b, 0 solve At x = b * * ICNTL(10) : number of iterative refinement when NRHS = 1 * * ICNTL(11) : > 0 return statistics * * ICNTL(12) : only used for SYM = 2, ordering strategy * * ICNTL(13) : * * ICNTL(14) : percentage of increase of the estimated working space * * ICNTL(15-17) : not used * * ICNTL(18) : only used if ICNTL(5) = 0, 0 matrix centralized, 1 structure on * host and mumps give the mapping, 2 structure on host and distributed matrix * for facto, 3 distributed matrix * * ICNTL(19) : > 0, Shur complement returned * * ICNTL(20) : 0 rhs dense, 1 rhs sparse * * ICNTL(21) : 0 solution in rhs, 1 solution distributed in ISOL_loc and SOL_loc * allocated by user * * ICNTL(22) : 0 in-core, 1 out-of-core * * ICNTL(23) : maximum memory allocatable by mumps pre proc * * ICNTL(24) : controls the detection of "null pivot rows" * * ICNTL(25) : * * ICNTL(26) : * * ICNTL(27) : * * ICNTL(28) : 0 automatic choice, 1 sequential analysis, 2 parallel analysis * * ICNTL(29) : 0 automatic choice, 1 PT-Scotch, 2 ParMetis */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "dof_manager_default.hh" #include "sparse_matrix_aij.hh" #if defined(AKANTU_USE_MPI) #include "static_communicator_mpi.hh" #include "mpi_type_wrapper.hh" #endif #include "solver_mumps.hh" #include "dof_synchronizer.hh" /* -------------------------------------------------------------------------- */ // static std::ostream & operator <<(std::ostream & stream, const DMUMPS_STRUC_C // & _this) { // stream << "DMUMPS Data [" << std::endl; // stream << " + job : " << _this.job << std::endl; // stream << " + par : " << _this.par << std::endl; // stream << " + sym : " << _this.sym << std::endl; // stream << " + comm_fortran : " << _this.comm_fortran << std::endl; // stream << " + nz : " << _this.nz << std::endl; // stream << " + irn : " << _this.irn << std::endl; // stream << " + jcn : " << _this.jcn << std::endl; // stream << " + nz_loc : " << _this.nz_loc << std::endl; // stream << " + irn_loc : " << _this.irn_loc << std::endl; // stream << " + jcn_loc : " << _this.jcn_loc << std::endl; // stream << "]"; // return stream; // } __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ SparseSolverMumps::SparseSolverMumps(DOFManagerDefault & dof_manager, const ID & matrix_id, const ID & id, const MemoryID & memory_id) : SparseSolver(dof_manager, matrix_id, id, memory_id), dof_manager(dof_manager), matrix(dof_manager.getMatrix(matrix_id)), rhs(dof_manager.getResidual()), solution(dof_manager.getGlobalSolution()), master_rhs_solution(0, 1) { AKANTU_DEBUG_IN(); StaticCommunicator & communicator = StaticCommunicator::getStaticCommunicator(); this->prank = communicator.whoAmI(); #ifdef AKANTU_USE_MPI this->parallel_method = _fully_distributed; #else // AKANTU_USE_MPI this->parallel_method = _not_parallel; #endif // AKANTU_USE_MPI this->mumps_data.par = 1; // The host is part of computations switch (this->parallel_method) { case _not_parallel: break; case _master_slave_distributed: this->mumps_data.par = 0; // The host is not part of the computations case _fully_distributed: #ifdef AKANTU_USE_MPI const StaticCommunicatorMPI & mpi_st_comm = dynamic_cast( communicator.getRealStaticCommunicator()); this->mumps_data.comm_fortran = MPI_Comm_c2f(mpi_st_comm.getMPITypeWrapper().getMPICommunicator()); #else AKANTU_DEBUG_ERROR( "You cannot use parallel method to solve without activating MPI"); #endif break; } this->mumps_data.sym = 2 * (this->matrix.getMatrixType() == _symmetric); this->prank = communicator.whoAmI(); this->mumps_data.job = _smj_initialize; // initialize dmumps_c(&this->mumps_data); /* ------------------------------------------------------------------------ */ this->setOutputLevel(); if (AKANTU_DEBUG_TEST(dblDump)) { strcpy(this->mumps_data.write_problem, "mumps_matrix.mtx"); } this->last_profile_release = this->matrix.getProfileRelease() - 1; this->last_value_release = this->matrix.getValueRelease() - 1; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SparseSolverMumps::~SparseSolverMumps() { AKANTU_DEBUG_IN(); this->mumps_data.job = _smj_destroy; // destroy dmumps_c(&this->mumps_data); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::setOutputLevel() { // Output setup icntl(1) = 0; // error output icntl(2) = 0; // dignostics output icntl(3) = 0; // informations icntl(4) = 0; #if !defined(AKANTU_NDEBUG) DebugLevel dbg_lvl = debug::debugger.getDebugLevel(); icntl(1) = (dbg_lvl >= dblWarning) ? 6 : 0; icntl(3) = (dbg_lvl >= dblInfo) ? 6 : 0; icntl(2) = (dbg_lvl >= dblTrace) ? 6 : 0; if (dbg_lvl >= dblDump) { icntl(4) = 4; } else if (dbg_lvl >= dblTrace) { icntl(4) = 3; } else if (dbg_lvl >= dblInfo) { icntl(4) = 2; } else if (dbg_lvl >= dblWarning) { icntl(4) = 1; } #endif } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::initMumpsData() { // Default Scaling icntl(8) = 77; // Assembled matrix icntl(5) = 0; /// Default centralized dense second member icntl(20) = 0; icntl(21) = 0; // automatic choice for analysis icntl(28) = 0; UInt size = this->matrix.getSize(); if (prank == 0) { this->master_rhs_solution.resize(size); } this->mumps_data.nz_alloc = 0; this->mumps_data.n = size; switch (this->parallel_method) { case _fully_distributed: icntl(18) = 3; // fully distributed this->mumps_data.nz_loc = matrix.getNbNonZero(); this->mumps_data.irn_loc = matrix.getIRN().storage(); this->mumps_data.jcn_loc = matrix.getJCN().storage(); break; case _not_parallel: case _master_slave_distributed: icntl(18) = 0; // centralized if (prank == 0) { this->mumps_data.nz = matrix.getNbNonZero(); this->mumps_data.irn = matrix.getIRN().storage(); this->mumps_data.jcn = matrix.getJCN().storage(); } else { this->mumps_data.nz = 0; this->mumps_data.irn = NULL; this->mumps_data.jcn = NULL; } break; default: AKANTU_DEBUG_ERROR("This case should not happen!!"); } } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::initialize() { AKANTU_DEBUG_IN(); this->analysis(); // icntl(14) = 80; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::analysis() { AKANTU_DEBUG_IN(); initMumpsData(); this->mumps_data.job = _smj_analyze; // analyze dmumps_c(&this->mumps_data); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::factorize() { AKANTU_DEBUG_IN(); this->mumps_data.rhs = this->rhs.storage(); if (parallel_method == _fully_distributed) this->mumps_data.a_loc = this->matrix.getA().storage(); else { if (prank == 0) this->mumps_data.a = this->matrix.getA().storage(); } this->mumps_data.job = _smj_factorize; // factorize dmumps_c(&this->mumps_data); this->printError(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::solve() { AKANTU_DEBUG_IN(); this->setOutputLevel(); this->dof_manager.updateGlobalBlockedDofs(); this->matrix.applyBoundary(); + this->master_rhs_solution.copy(this->dof_manager.getResidual()); // if (prank == 0) { // matrix.getDOFSynchronizer().gather(this->rhs, 0, this->master_rhs_solution); // } else { // this->matrix.getDOFSynchronizer().gather(this->rhs, 0); // } if (this->last_profile_release != this->matrix.getProfileRelease()) { this->analysis(); this->last_profile_release = this->matrix.getProfileRelease(); } if (this->last_value_release != this->matrix.getValueRelease()) { this->factorize(); this->last_value_release = this->matrix.getValueRelease(); } - if (prank == 0) + if (prank == 0) { this->mumps_data.rhs = this->master_rhs_solution.storage(); + } this->mumps_data.job = _smj_solve; // solve dmumps_c(&this->mumps_data); this->printError(); + this->dof_manager.getGlobalSolution().copy(this->master_rhs_solution); // if (prank == 0) { // matrix.getDOFSynchronizer().scatter(this->solution, 0, this->master_rhs_solution); // } else { // this->matrix.getDOFSynchronizer().gather(this->solution, 0); // } + this->dof_manager.splitSolutionPerDOFs(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseSolverMumps::printError() { Vector _info_v(2); _info_v[0] = info(1); // to get errors _info_v[1] = -info(1); // to get warnings StaticCommunicator::getStaticCommunicator().allReduce(_info_v, _so_min); _info_v[1] = -_info_v[1]; if (_info_v[0] < 0) { // < 0 is an error switch (_info_v[0]) { case -10: AKANTU_DEBUG_ERROR("The matrix is singular"); break; case -9: { icntl(14) += 10; if (icntl(14) != 90) { // std::cout << "Dynamic memory increase of 10%" << std::endl; AKANTU_DEBUG_WARNING("MUMPS dynamic memory is insufficient it will be " "increased allowed to use 10% more"); // change releases to force a recompute this->last_value_release--; this->last_profile_release--; this->solve(); } else { AKANTU_DEBUG_ERROR("The MUMPS workarray is too small INFO(2)=" << info(2) << "No further increase possible"); break; } } default: AKANTU_DEBUG_ERROR("Error in mumps during solve process, check mumps " "user guide INFO(1) = " << _info_v[1]); } } else if (_info_v[1] > 0) { AKANTU_DEBUG_WARNING("Warning in mumps during solve process, check mumps " "user guide INFO(1) = " << _info_v[1]); } } __END_AKANTU__ diff --git a/src/solver/sparse_matrix.cc b/src/solver/sparse_matrix.cc index 055c2929e..4b9069d7c 100644 --- a/src/solver/sparse_matrix.cc +++ b/src/solver/sparse_matrix.cc @@ -1,74 +1,82 @@ /** * @file sparse_matrix.cc * * @author Nicolas Richart * * @date creation: Mon Dec 13 2010 * @date last modification: Mon Nov 16 2015 * * @brief implementation of the SparseMatrix class * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include "sparse_matrix.hh" #include "static_communicator.hh" #include "dof_manager.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ SparseMatrix::SparseMatrix(DOFManager & dof_manager, const MatrixType & matrix_type, const ID & id) : id(id), _dof_manager(dof_manager), matrix_type(matrix_type), size(dof_manager.getSystemSize()), nb_non_zero(0), matrix_release(1) { AKANTU_DEBUG_IN(); StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); this->nb_proc = comm.getNbProc(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SparseMatrix::SparseMatrix(const SparseMatrix & matrix, const ID & id) : id(id), _dof_manager(matrix._dof_manager), matrix_type(matrix.matrix_type), size(matrix.size), nb_proc(matrix.nb_proc), nb_non_zero(matrix.nb_non_zero), matrix_release(1) {} /* -------------------------------------------------------------------------- */ SparseMatrix::~SparseMatrix() {} /* -------------------------------------------------------------------------- */ Array & operator*=(Array & vect, const SparseMatrix & mat) { - Array tmp(vect.getSize(), vect.getNbComponent()); + Array tmp(vect.getSize(), vect.getNbComponent(), 0.); mat.matVecMul(vect, tmp); vect.copy(tmp); return vect; } +/* -------------------------------------------------------------------------- */ +void SparseMatrix::add(const SparseMatrix & B, Real alpha) { + B.addMeTo(*this, alpha); +} + +/* -------------------------------------------------------------------------- */ + + __END_AKANTU__ diff --git a/src/solver/sparse_matrix.hh b/src/solver/sparse_matrix.hh index 5ac8300d0..de6adfe4c 100644 --- a/src/solver/sparse_matrix.hh +++ b/src/solver/sparse_matrix.hh @@ -1,148 +1,158 @@ /** * @file sparse_matrix.hh * * @author Nicolas Richart * * @date creation: Mon Dec 13 2010 * @date last modification: Fri Oct 16 2015 * * @brief sparse matrix storage class (distributed assembled matrix) * This is a COO format (Coordinate List) * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SPARSE_MATRIX_HH__ #define __AKANTU_SPARSE_MATRIX_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class DOFManager; } __BEGIN_AKANTU__ class SparseMatrix { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: SparseMatrix(DOFManager & dof_manager, const MatrixType & matrix_type, const ID & id = "sparse_matrix"); SparseMatrix(const SparseMatrix & matrix, const ID & id = "sparse_matrix"); virtual ~SparseMatrix(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// remove the existing profile virtual void clearProfile(); /// set the matrix to 0 virtual void clear() = 0; /// add a non-zero element to the profile virtual inline UInt addToProfile(UInt i, UInt j) = 0; /// assemble a local matrix in the sparse one - virtual inline void addToMatrix(UInt i, UInt j, Real value) = 0; + virtual void addToMatrix(UInt i, UInt j, Real value) = 0; /// save the profil in a file using the MatrixMarket file format virtual void saveProfile(const std::string &) const { AKANTU_DEBUG_TO_IMPLEMENT(); } /// save the matrix in a file using the MatrixMarket file format virtual void saveMatrix(const std::string &) const { AKANTU_DEBUG_TO_IMPLEMENT(); }; - /// add matrix assuming the profile are the same - virtual void add(const SparseMatrix & matrix, Real alpha = 1.) = 0; + /// multiply the matrix by a coefficient + virtual void mul(Real alpha) = 0; + + /// add matrices + virtual void add(const SparseMatrix & matrix, Real alpha = 1.); /// Equivalent of *gemv in blas virtual void matVecMul(const Array & x, Array & y, Real alpha = 1., Real beta = 0.) const = 0; /// modify the matrix to "remove" the blocked dof virtual void applyBoundary(Real block_val = 1.) = 0; + /// operator *= + SparseMatrix & operator*=(Real alpha) { this->mul(alpha); return *this; } + +protected: + /// This is the revert of add B += \alpha * *this; + virtual void addMeTo(SparseMatrix & B, Real alpha) const = 0; + /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the values at potition i, j virtual inline Real operator()(UInt i, UInt j) const { AKANTU_DEBUG_TO_IMPLEMENT(); } /// return the values at potition i, j virtual inline Real & operator()(UInt i, UInt j) { AKANTU_DEBUG_TO_IMPLEMENT(); } AKANTU_GET_MACRO(NbNonZero, nb_non_zero, UInt); AKANTU_GET_MACRO(Size, size, UInt); AKANTU_GET_MACRO(MatrixType, matrix_type, const MatrixType &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: ID id; /// Underlying dof manager DOFManager & _dof_manager; /// sparce matrix type MatrixType matrix_type; /// Size of the matrix UInt size; /// number of processors UInt nb_proc; /// number of non zero element UInt nb_non_zero; /// matrix definition releasez UInt matrix_release; }; Array & operator*=(Array & vect, const SparseMatrix & mat); __END_AKANTU__ /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "sparse_matrix_inline_impl.cc" #endif /* __AKANTU_SPARSE_MATRIX_HH__ */ diff --git a/src/solver/sparse_matrix_aij.cc b/src/solver/sparse_matrix_aij.cc index 32c5a8d9e..cf113e09a 100644 --- a/src/solver/sparse_matrix_aij.cc +++ b/src/solver/sparse_matrix_aij.cc @@ -1,308 +1,225 @@ /** * @file sparse_matrix_aij.cc * * @author Nicolas Richart * * @date Tue Aug 18 16:31:23 2015 * * @brief Implementation of the AIJ sparse matrix * * @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 "sparse_matrix_aij.hh" #include "dof_manager_default.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ SparseMatrixAIJ::SparseMatrixAIJ(DOFManagerDefault & dof_manager, const MatrixType & matrix_type, const ID & id) : SparseMatrix(dof_manager, matrix_type, id), dof_manager(dof_manager), irn(0, 1, id + ":irn"), jcn(0, 1, id + ":jcn"), a(0, 1, id + ":a"), profile_release(0), value_release(0) {} /* -------------------------------------------------------------------------- */ SparseMatrixAIJ::SparseMatrixAIJ(const SparseMatrixAIJ & matrix, const ID & id) : SparseMatrix(matrix, id), dof_manager(matrix.dof_manager), irn(matrix.irn, true, id + ":irn"), jcn(matrix.jcn, true, id + ":jcn"), a(matrix.a, true, id + ":a"), profile_release(0), value_release(0) {} /* -------------------------------------------------------------------------- */ SparseMatrixAIJ::~SparseMatrixAIJ() {} -// /* -------------------------------------------------------------------------- -// */ -// void SparseMatrixAIJ::buildProfile(const Mesh & mesh, -// const DOFSynchronizer & dof_synchronizer, -// UInt nb_degree_of_freedom) { -// AKANTU_DEBUG_IN(); - -// // if(irn_jcn_to_k) delete irn_jcn_to_k; -// // irn_jcn_to_k = new std::map, UInt>; -// clearProfile(); - -// this->dof_synchronizer = &const_cast(dof_synchronizer); - -// coordinate_list_map::iterator irn_jcn_k_it; - -// Int * eq_nb_val = dof_synchronizer.getGlobalDOFEquationNumbers().storage(); - -// Mesh::type_iterator it = -// mesh.firstType(mesh.getSpatialDimension(), _not_ghost, -// _ek_not_defined); -// Mesh::type_iterator end = -// mesh.lastType(mesh.getSpatialDimension(), _not_ghost, _ek_not_defined); -// for (; it != end; ++it) { -// UInt nb_element = mesh.getNbElement(*it); -// UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); -// UInt size_mat = nb_nodes_per_element * nb_degree_of_freedom; - -// UInt * conn_val = mesh.getConnectivity(*it, _not_ghost).storage(); -// Int * local_eq_nb_val = -// new Int[nb_degree_of_freedom * nb_nodes_per_element]; - -// for (UInt e = 0; e < nb_element; ++e) { -// Int * tmp_local_eq_nb_val = local_eq_nb_val; -// for (UInt i = 0; i < nb_nodes_per_element; ++i) { -// UInt n = conn_val[i]; -// for (UInt d = 0; d < nb_degree_of_freedom; ++d) { -// *tmp_local_eq_nb_val++ = eq_nb_val[n * nb_degree_of_freedom + d]; -// } -// // memcpy(tmp_local_eq_nb_val, eq_nb_val + n * nb_degree_of_freedom, -// // nb_degree_of_freedom * sizeof(Int)); -// // tmp_local_eq_nb_val += nb_degree_of_freedom; -// } - -// for (UInt i = 0; i < size_mat; ++i) { -// UInt c_irn = local_eq_nb_val[i]; -// if (c_irn < this->size) { -// UInt j_start = (sparse_matrix_type == _symmetric) ? i : 0; -// for (UInt j = j_start; j < size_mat; ++j) { -// UInt c_jcn = local_eq_nb_val[j]; -// if (c_jcn < this->size) { -// KeyCOO irn_jcn = key(c_irn, c_jcn); -// irn_jcn_k_it = irn_jcn_k.find(irn_jcn); - -// if (irn_jcn_k_it == irn_jcn_k.end()) { -// irn_jcn_k[irn_jcn] = nb_non_zero; -// irn.push_back(c_irn + 1); -// jcn.push_back(c_jcn + 1); -// nb_non_zero++; -// } -// } -// } -// } -// } -// conn_val += nb_nodes_per_element; -// } - -// delete[] local_eq_nb_val; -// } - -// /// for pbc @todo correct it for parallel -// if (StaticCommunicator::getStaticCommunicator().getNbProc() == 1) { -// for (UInt i = 0; i < size; ++i) { -// KeyCOO irn_jcn = key(i, i); -// irn_jcn_k_it = irn_jcn_k.find(irn_jcn); -// if (irn_jcn_k_it == irn_jcn_k.end()) { -// irn_jcn_k[irn_jcn] = nb_non_zero; -// irn.push_back(i + 1); -// jcn.push_back(i + 1); -// nb_non_zero++; -// } -// } -// } - -// a.resize(nb_non_zero); - -// AKANTU_DEBUG_OUT(); -// } - /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::applyBoundary(Real block_val) { AKANTU_DEBUG_IN(); const Array & blocked_dofs = this->dof_manager.getGlobalBlockedDOFs(); Array::const_scalar_iterator irn_val = irn.begin(); Array::const_scalar_iterator jcn_val = jcn.begin(); Array::scalar_iterator a_val = a.begin(); for (UInt i = 0; i < nb_non_zero; ++i) { UInt ni = this->dof_manager.globalToLocalEquationNumber(*irn_val - 1); UInt nj = this->dof_manager.globalToLocalEquationNumber(*jcn_val - 1); if (blocked_dofs(ni) || blocked_dofs(nj)) { if (*irn_val != *jcn_val) { *a_val = 0; } else { if (this->dof_manager.isLocalOrMasterDOF(ni)) { *a_val = block_val; } else { *a_val = 0; } } } ++irn_val; ++jcn_val; ++a_val; } this->value_release++; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::saveProfile(const std::string & filename) const { AKANTU_DEBUG_IN(); std::ofstream outfile; outfile.open(filename.c_str()); outfile << "%%MatrixMarket matrix coordinate pattern"; if (this->matrix_type == _symmetric) outfile << " symmetric"; else outfile << " general"; outfile << std::endl; UInt m = this->size; outfile << m << " " << m << " " << this->nb_non_zero << std::endl; for (UInt i = 0; i < this->nb_non_zero; ++i) { outfile << this->irn.storage()[i] << " " << this->jcn.storage()[i] << " 1" << std::endl; } outfile.close(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::saveMatrix(const std::string & filename) const { AKANTU_DEBUG_IN(); std::ofstream outfile; outfile.precision(std::numeric_limits::digits10); outfile.open(filename.c_str()); outfile << "%%MatrixMarket matrix coordinate real"; if (this->matrix_type == _symmetric) outfile << " symmetric"; else outfile << " general"; outfile << std::endl; outfile << this->size << " " << this->size << " " << this->nb_non_zero << std::endl; for (UInt i = 0; i < this->nb_non_zero; ++i) { outfile << this->irn(i) << " " << this->jcn(i) << " " << this->a(i) << std::endl; } outfile.close(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::matVecMul(const Array & x, Array & y, Real alpha, Real beta) const { AKANTU_DEBUG_IN(); y *= beta; Array::const_scalar_iterator i_it = this->irn.storage(); Array::const_scalar_iterator j_it = this->jcn.storage(); Array::const_scalar_iterator a_it = this->a.storage(); Array::const_scalar_iterator x_it = x.storage(); Array::scalar_iterator y_it = y.storage(); for (UInt k = 0; k < this->nb_non_zero; ++k, ++i_it, ++j_it, ++a_it) { - const Int & i = *i_it; - const Int & j = *j_it; + Int i = *i_it - 1; + Int j = *j_it - 1; const Real & A = *a_it; y_it[i] += alpha * A * x_it[j]; if ((this->matrix_type == _symmetric) && (i != j)) y_it[j] += alpha * A * x_it[i]; } // if (dof_synchronizer) // dof_synchronizer->reduceSynchronize(vect); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::copyContent(const SparseMatrix & matrix) { AKANTU_DEBUG_IN(); const SparseMatrixAIJ & mat = dynamic_cast(matrix); AKANTU_DEBUG_ASSERT(nb_non_zero == mat.getNbNonZero(), "The to matrix don't have the same profiles"); memcpy(a.storage(), mat.getA().storage(), nb_non_zero * sizeof(Real)); this->value_release++; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -void SparseMatrixAIJ::add(const SparseMatrix & B, Real alpha) { - Array::scalar_iterator a_it = this->a.begin(); - - try { - const SparseMatrixAIJ & B_aij = dynamic_cast(B); - Array::const_scalar_iterator b_it = B_aij.a.begin(); - for (UInt n = 0; n < this->nb_non_zero; ++n, ++a_it, ++b_it) { - *a_it += alpha * *b_it; - } - } catch(...) { - Array::const_scalar_iterator i_it = this->irn.begin(); - Array::const_scalar_iterator j_it = this->jcn.begin(); - for (UInt n = 0; n < this->nb_non_zero; ++n, ++a_it, ++i_it, ++j_it) { - const Int & i = *i_it; - const Int & j = *j_it; - Real & A_ij = *a_it; - A_ij += alpha * B(i - 1, j - 1); - } +template +void SparseMatrixAIJ::addMeToTemplated(MatrixType & B, Real alpha) const { + Array::const_scalar_iterator i_it = this->irn.begin(); + Array::const_scalar_iterator j_it = this->jcn.begin(); + Array::const_scalar_iterator a_it = this->a.begin(); + for (UInt n = 0; n < this->nb_non_zero; ++n, ++a_it, ++i_it, ++j_it) { + const Int & i = *i_it; + const Int & j = *j_it; + const Real & A_ij = *a_it; + B.addToMatrix(i - 1, j - 1, alpha * A_ij); } +} + +/* -------------------------------------------------------------------------- */ +void SparseMatrixAIJ::addMeTo(SparseMatrix & B, Real alpha) const { + if(SparseMatrixAIJ * B_aij = dynamic_cast(&B)) { + this->addMeToTemplated(*B_aij, alpha); + } else { + // this->addMeToTemplated(*this, alpha); + } +} +/* -------------------------------------------------------------------------- */ +void SparseMatrixAIJ::mul(Real alpha) { + this->a *= alpha; this->value_release++; } /* -------------------------------------------------------------------------- */ void SparseMatrixAIJ::clear() { a.set(0.); this->value_release++; } __END_AKANTU__ diff --git a/src/solver/sparse_matrix_aij.hh b/src/solver/sparse_matrix_aij.hh index 2c32b083b..79c07d4bf 100644 --- a/src/solver/sparse_matrix_aij.hh +++ b/src/solver/sparse_matrix_aij.hh @@ -1,165 +1,178 @@ /** * @file sparse_matrix_aij.hh * * @author Nicolas Richart * * @date Mon Aug 17 21:22:57 2015 * * @brief AIJ implementation of the SparseMatrix (this the format used by Mumps) * * @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 "sparse_matrix.hh" +#include "aka_common.hh" #include "aka_array.hh" +#include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SPARSE_MATRIX_AIJ_HH__ #define __AKANTU_SPARSE_MATRIX_AIJ_HH__ namespace akantu { -class DOFManagerDefault; + class DOFManagerDefault; } __BEGIN_AKANTU__ class SparseMatrixAIJ : public SparseMatrix { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: SparseMatrixAIJ(DOFManagerDefault & dof_manager, const MatrixType & matrix_type, const ID & id = "sparse_matrix_aij"); SparseMatrixAIJ(const SparseMatrixAIJ & matrix, const ID & id = "sparse_matrix_aij"); virtual ~SparseMatrixAIJ(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// remove the existing profile inline void clearProfile(); /// add a non-zero element virtual UInt addToProfile(UInt i, UInt j); /// set the matrix to 0 virtual void clear(); /// assemble a local matrix in the sparse one inline void addToMatrix(UInt i, UInt j, Real value); /// set the size of the matrix void resize(UInt size) { this->size = size; } /// modify the matrix to "remove" the blocked dof virtual void applyBoundary(Real block_val = 1.); /// save the profil in a file using the MatrixMarket file format virtual void saveProfile(const std::string & filename) const; /// save the matrix in a file using the MatrixMarket file format virtual void saveMatrix(const std::string & filename) const; /// copy assuming the profile are the same virtual void copyContent(const SparseMatrix & matrix); - /// add matrix assuming the profile are the same - virtual void add(const SparseMatrix & matrix, Real alpha); + /// multiply the matrix by a scalar + void mul(Real alpha); + + /// add matrix *this += B + //virtual void add(const SparseMatrix & matrix, Real alpha); /// Equivalent of *gemv in blas virtual void matVecMul(const Array & x, Array & y, Real alpha = 1., Real beta = 0.) const; /* ------------------------------------------------------------------------ */ /// accessor to A_{ij} - if (i, j) not present it returns 0 inline Real operator()(UInt i, UInt j) const; /// accessor to A_{ij} - if (i, j) not present it fails, (i, j) should be /// first added to the profile inline Real & operator()(UInt i, UInt j); +protected: + /// This is the revert of add B += \alpha * *this; + virtual void addMeTo(SparseMatrix & B, Real alpha) const; + +private: + /// This is just to inline the addToMatrix function + template + void addMeToTemplated(MatrixType & B, Real alpha) const; + /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(IRN, irn, const Array &); AKANTU_GET_MACRO(JCN, jcn, const Array &); AKANTU_GET_MACRO(A, a, const Array &); /// The release changes at each call of a function that changes the profile, /// it in increasing but could overflow so it should be checked as /// (my_release != release) and not as (my_release < release) AKANTU_GET_MACRO(ProfileRelease, profile_release, UInt); AKANTU_GET_MACRO(ValueRelease, value_release, UInt); protected: typedef std::pair KeyCOO; typedef unordered_map::type coordinate_list_map; /// get the pair corresponding to (i, j) inline KeyCOO key(UInt i, UInt j) const { if (this->matrix_type == _symmetric && (i > j)) return std::make_pair(j, i); return std::make_pair(i, j); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: DOFManagerDefault & dof_manager; /// row indexes Array irn; /// column indexes Array jcn; /// values : A[k] = Matrix[irn[k]][jcn[k]] Array a; /// Profile release UInt profile_release; /// Value release UInt value_release; /* * map for (i, j) -> k correspondence */ coordinate_list_map irn_jcn_k; }; __END_AKANTU__ /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "sparse_matrix_aij_inline_impl.cc" #endif /* __AKANTU_SPARSE_MATRIX_AIJ_HH__ */ diff --git a/src/solver/sparse_matrix_aij_inline_impl.cc b/src/solver/sparse_matrix_aij_inline_impl.cc index 5a084fafd..17d4ded0d 100644 --- a/src/solver/sparse_matrix_aij_inline_impl.cc +++ b/src/solver/sparse_matrix_aij_inline_impl.cc @@ -1,115 +1,111 @@ /** * @file sparse_matrix_aij_inline_impl.cc * * @author Nicolas Richart * * @date Tue Aug 18 00:42:45 2015 * * @brief Implementation of inline functions of SparseMatrixAIJ * * @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_SPARSE_MATRIX_AIJ_INLINE_IMPL_CC__ #define __AKANTU_SPARSE_MATRIX_AIJ_INLINE_IMPL_CC__ __BEGIN_AKANTU__ inline UInt SparseMatrixAIJ::addToProfile(UInt i, UInt j) { KeyCOO jcn_irn = this->key(i, j); coordinate_list_map::iterator it = this->irn_jcn_k.find(jcn_irn); + if (!(it == this->irn_jcn_k.end())) return it->second; if(i + 1 > this->size) this->size = i + 1; if(j + 1 > this->size) this->size = j + 1; this->irn.push_back(i + 1); this->jcn.push_back(j + 1); this->a.push_back(0.); this->irn_jcn_k[jcn_irn] = this->nb_non_zero; (this->nb_non_zero)++; this->profile_release++; this->value_release++; return (this->nb_non_zero - 1); } /* -------------------------------------------------------------------------- */ inline void SparseMatrixAIJ::clearProfile() { SparseMatrix::clearProfile(); this->irn_jcn_k.clear(); this->irn.resize(0); this->jcn.resize(0); this->a.resize(0); this->size = 0; this->profile_release++; this->value_release++; } /* -------------------------------------------------------------------------- */ inline void SparseMatrixAIJ::addToMatrix(UInt i, UInt j, Real value) { - KeyCOO jcn_irn = this->key(i, j); - coordinate_list_map::iterator irn_jcn_k_it = this->irn_jcn_k.find(jcn_irn); - - AKANTU_DEBUG_ASSERT(irn_jcn_k_it != this->irn_jcn_k.end(), - "Couple (i,j) = (" << i << "," << j - << ") does not exist in the profile"); + UInt idx = this->addToProfile(i, j); - this->a(irn_jcn_k_it->second) += value; + this->a(idx) += value; this->value_release++; } /* -------------------------------------------------------------------------- */ inline Real SparseMatrixAIJ::operator()(UInt i, UInt j) const { KeyCOO jcn_irn = this->key(i, j); coordinate_list_map::const_iterator irn_jcn_k_it = this->irn_jcn_k.find(jcn_irn); if (irn_jcn_k_it == this->irn_jcn_k.end()) return 0.; return this->a(irn_jcn_k_it->second); } /* -------------------------------------------------------------------------- */ inline Real & SparseMatrixAIJ::operator()(UInt i, UInt j) { KeyCOO jcn_irn = this->key(i, j); coordinate_list_map::iterator irn_jcn_k_it = this->irn_jcn_k.find(jcn_irn); AKANTU_DEBUG_ASSERT(irn_jcn_k_it != this->irn_jcn_k.end(), "Couple (i,j) = (" << i << "," << j << ") does not exist in the profile"); // it may change the profile so it is considered as a change this->value_release++; return this->a(irn_jcn_k_it->second); } __END_AKANTU__ #endif /* __AKANTU_SPARSE_MATRIX_AIJ_INLINE_IMPL_CC__ */ diff --git a/src/synchronizer/dof_synchronizer.hh b/src/synchronizer/dof_synchronizer.hh index 90acb22ef..ecea7691b 100644 --- a/src/synchronizer/dof_synchronizer.hh +++ b/src/synchronizer/dof_synchronizer.hh @@ -1,245 +1,244 @@ /** * @file dof_synchronizer.hh * * @author Aurelia Isabel Cuba Ramos * @author Nicolas Richart * * @date creation: Fri Jun 17 2011 * @date last modification: Tue Dec 08 2015 * * @brief Synchronize Array of DOFs * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_array.hh" #include "static_communicator.hh" #include "synchronizer.hh" /* -------------------------------------------------------------------------- */ - - +namespace akantu { + class Mesh; +} #ifndef __AKANTU_DOF_SYNCHRONIZER_HH__ #define __AKANTU_DOF_SYNCHRONIZER_HH__ __BEGIN_AKANTU__ -class Mesh; - template class AddOperation { public: inline T operator()(T & b, T & a) { return a + b; }; }; class DOFSynchronizer : public Synchronizer { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: DOFSynchronizer(const Mesh & mesh, UInt nb_degree_of_freedom); virtual ~DOFSynchronizer(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// asynchronous synchronization of ghosts virtual void asynchronousSynchronize(DataAccessor & data_accessor,SynchronizationTag tag); /// wait end of asynchronous synchronization of ghosts virtual void waitEndSynchronize(DataAccessor & data_accessor,SynchronizationTag tag); /// compute buffer size for a given tag and data accessor virtual void computeBufferSize(DataAccessor & data_accessor, SynchronizationTag tag); /// init the scheme for scatter and gather operation, need extra memory void initScatterGatherCommunicationScheme(); /// initialize the equation number with local ids void initLocalDOFEquationNumbers(); /// initialize the equation number with local ids void initGlobalDOFEquationNumbers(); /** * Gather the DOF value on the root proccessor * * @param to_gather data to gather * @param root processor on which data are gathered * @param gathered Array containing the gathered data, only valid on root processor */ template /// Gather the DOF value on the root proccessor void gather(const Array & to_gather, UInt root, Array * gathered = NULL) const; /** * Scatter a DOF Array form root to all processors * * @param scattered data to scatter, only valid on root processor * @param root processor scattering data * @param to_scatter result of scattered data */ template /// Scatter a DOF Array form root to all processors void scatter(Array & scattered, UInt root, const Array * to_scatter = NULL) const; template void synchronize(Array & vector) const ; template