diff --git a/CMakeLists.txt b/CMakeLists.txt index f476be652..8c723f373 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,190 +1,194 @@ #=============================================================================== # @file CMakeLists.txt # # @author Guillaume Anciaux # @author Nicolas Richart # # @date creation: Mon Jun 14 2010 # @date last modification: Fri Jan 22 2016 # # @brief main configuration file # # @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 #------------------------------------------------------------------------------- # _ _ # | | | | # __ _| | ____ _ _ __ | |_ _ _ # / _` | |/ / _` | '_ \| __| | | | # | (_| | < (_| | | | | |_| |_| | # \__,_|_|\_\__,_|_| |_|\__|\__,_| # #=============================================================================== #=============================================================================== # CMake Project #=============================================================================== cmake_minimum_required(VERSION 3.1.3) # add this options before PROJECT keyword set(CMAKE_DISABLE_SOURCE_CHANGES ON) set(CMAKE_DISABLE_IN_SOURCE_BUILD ON) project(Akantu) enable_language(CXX) #=============================================================================== # Misc. config for cmake #=============================================================================== set(AKANTU_CMAKE_DIR "${PROJECT_SOURCE_DIR}/cmake") set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/cmake") set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/cmake/Modules") set(BUILD_SHARED_LIBS ON CACHE BOOL "Build shared libraries.") mark_as_advanced(BUILD_SHARED_LIBS) if(NOT AKANTU_TARGETS_EXPORT) set(AKANTU_TARGETS_EXPORT AkantuTargets) endif() include(CMakeVersionGenerator) include(CMakePackagesSystem) include(CMakeFlagsHandling) include(AkantuPackagesSystem) include(AkantuMacros) include(AkantuCleaning) #cmake_activate_debug_message() #=============================================================================== # Version Number #=============================================================================== # AKANTU version number. An even minor number corresponds to releases. set(AKANTU_MAJOR_VERSION 3) set(AKANTU_MINOR_VERSION 0) set(AKANTU_PATCH_VERSION 0) define_project_version() #=============================================================================== # Options #=============================================================================== # Debug set(CMAKE_CXX_FLAGS_RELEASE "-O3 -DNDEBUG -DAKANTU_NDEBUG" CACHE STRING "Flags used by the compiler during release builds" FORCE) #add_flags(cxx "-Wall -Wextra -pedantic -Werror") if (CMAKE_CXX_COMPILER_ID STREQUAL "GNU" OR CMAKE_CXX_COMPILER_ID STREQUAL "Clang") add_flags(cxx "-Wall -Wextra -pedantic") # -Weffc++ + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG_INIT} -ggdb3" + CACHE STRING "Flags used by the compiler during debug builds" FORCE) + set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "${CMAKE_CXX_FLAGS_RELWITHDEBINFO_INIT} -ggdb3" + CACHE STRING "Flags used by the compiler during debug builds" FORCE) else() add_flags(cxx "-Wall") endif() option(AKANTU_EXAMPLES "Activate examples" OFF) option(AKANTU_TESTS "Activate tests" OFF) set(AKANTU_PREFERRED_PYTHON_VERSION 3 CACHE STRING "Preferred version for python related things") mark_as_advanced(AKANTU_PREFERRED_PYTHON_VERSION) include(AkantuExtraCompilationProfiles) #=============================================================================== # Dependencies #=============================================================================== declare_akantu_types() package_list_packages(${PROJECT_SOURCE_DIR}/packages EXTRA_PACKAGES_FOLDER ${PROJECT_SOURCE_DIR}/extra_packages NO_AUTO_COMPILE_FLAGS) ## meta option \todo better way to do it when multiple package give enable the ## same feature if(AKANTU_SCOTCH) set(AKANTU_PARTITIONER ON) else() set(AKANTU_PARTITIONER OFF) endif() if(AKANTU_MUMPS) set(AKANTU_SOLVER ON) else() set(AKANTU_SOLVER OFF) endif() #=============================================================================== # Akantu library #=============================================================================== add_subdirectory(src) #=============================================================================== # Documentation #=============================================================================== if(AKANTU_DOCUMENTATION_DOXYGEN OR AKANTU_DOCUMENTATION_MANUAL) add_subdirectory(doc) else() set(AKANTU_DOC_EXCLUDE_FILES "${PROJECT_SOURCE_DIR}/doc/manual" CACHE INTERNAL "") endif() #=============================================================================== # Examples and tests #=============================================================================== include(AkantuTestsMacros) include(AkantuExampleMacros) if(AKANTU_TESTS) option(AKANTU_BUILD_ALL_TESTS "Build all tests" ON) find_package(GMSH REQUIRED) # package_is_activated(pybind11 _pybind11_act) # if(_pybind11_act) # find_package(pybind11 CONFIG REQUIRED QUIET) # to get the pybind11_add_module macro # endif() endif() if(AKANTU_EXAMPLES) find_package(GMSH REQUIRED) add_subdirectory(examples) endif() # tests add_test_tree(test) #=============================================================================== # Python interface #=============================================================================== package_is_activated(python_interface _python_act) if(_python_act) if(IS_ABSOLUTE "${CMAKE_INSTALL_PREFIX}") set(AKANTU_PYTHON_INSTALL_PREFIX ${CMAKE_INSTALL_PREFIX}) else() set(AKANTU_PYTHON_INSTALL_PREFIX "${PROJECT_BINARY_DIR}/${CMAKE_INSTALL_PREFIX}") endif() add_subdirectory(python) endif() #=============================================================================== # Install and Packaging #=============================================================================== include(AkantuInstall) option(AKANTU_DISABLE_CPACK "This option commands the generation of extra info for the \"make package\" target" ON) mark_as_advanced(AKANTU_DISABLE_CPACK) if(NOT AKANTU_DISABLE_CPACK) include(AkantuCPack) endif() diff --git a/Jenkinsfile b/Jenkinsfile index f33b6f40d..cfa30bfd5 100644 --- a/Jenkinsfile +++ b/Jenkinsfile @@ -1,78 +1,128 @@ pipeline { + parameters {string(defaultValue: '', description: 'api-token', name: 'API_TOKEN') + string(defaultValue: '', description: 'buildable phid', name: 'TARGET_PHID') + string(defaultValue: '', description: 'Commit id', name: 'COMMIT_ID') + string(defaultValue: '', description: 'Diff id', name: 'DIFF_ID') + } + + options { + disableConcurrentBuilds() + } + agent { dockerfile { additionalBuildArgs '--tag akantu-environment' } } stages { stage('Configure') { steps { - sh 'env' - sh 'mkdir -p build' - sh 'cd build; cmake -DAKANTU_COHESIVE_ELEMENT:BOOL=TRUE -DAKANTU_IMPLICIT:BOOL=TRUE -DAKANTU_PARALLEL:BOOL=TRUE -DAKANTU_PYTHON_INTERFACE:BOOL=TRUE -DAKANTU_TESTS:BOOL=TRUE ..' + sh """ + env + mkdir -p build + cd build + cmake -DAKANTU_COHESIVE_ELEMENT:BOOL=TRUE -DAKANTU_IMPLICIT:BOOL=TRUE -DAKANTU_PARALLEL:BOOL=TRUE -DAKANTU_PYTHON_INTERFACE:BOOL=TRUE -DAKANTU_TESTS:BOOL=TRUE .. + """ } } stage('Compile') { steps { sh 'make -C build/src || true' } } stage ('Warnings gcc') { steps { warnings(consoleParsers: [[parserName: 'GNU Make + GNU C Compiler (gcc)']]) } } stage('Compile python') { steps { sh 'make -C build/python || true' } } stage('Compile tests') { steps { sh 'make -C build/test || true' } } stage('Tests') { steps { - sh 'rm -rf build/gtest_reports' - sh 'cd build/ && ctest -T test --no-compress-output || true' - sh 'cp build/Testing/`head -n 1 < build/Testing/TAG`/Test.xml CTestResults.xml' + sh """ + rm -rf build/gtest_reports + cd build/ + #source ./akantu_environement.sh + ctest -T test --no-compress-output || true + TAG=`head -n 1 < build/Testing/TAG` + if [ -e build/Testing/${TAG}/Test.xml ]; then + cp build/Testing/${TAG}/Test.xml CTestResults.xml + fi + """ } } } environment { BLA_VENDOR = 'OpenBLAS' OMPI_MCA_plm = 'isolated' OMPI_MCA_btl = 'tcp,self' } post { always { step([$class: 'XUnitBuilder', thresholds: [ [$class: 'SkippedThreshold', failureThreshold: '0'], [$class: 'FailedThreshold', failureThreshold: '0']], tools: [[$class: 'CTestType', pattern: 'CTestResults.xml']]]) step([$class: 'XUnitBuilder', thresholds: [ [$class: 'SkippedThreshold', failureThreshold: '100'], [$class: 'FailedThreshold', failureThreshold: '0']], - tools: [[$class: 'GoogleTestType', pattern: 'build/gtest_reports/**']]]) + tools: [[$class: 'GoogleTestType', pattern: 'build/gtest_reports/**']]]) + archiveArtifacts artifacts: 'build/Testing/**', fingerprint: true + createartifact() + } + + success { + send_fail_pass('pass') } failure { emailext( body: '''${SCRIPT, template="groovy-html.template"}''', mimeType: 'text/html', subject: "[Jenkins] ${currentBuild.fullDisplayName} Failed", recipientProviders: [[$class: 'CulpritsRecipientProvider']], to: 'akantu-admins@akantu.ch', replyTo: 'akantu-admins@akantu.ch', attachLog: true, compressLog: false) + send_fail_pass('fail') } } } + +def send_fail_pass(state) { + sh """ +set +x +curl https://c4science.ch/api/harbormaster.sendmessage \ +-d api.token=${API_TOKEN} \ +-d buildTargetPHID=${TARGET_PHID} \ +-d type=${state} +""" +} + +def createartifact() { + sh """ set +x +curl https://c4science.ch/api/harbormaster.createartifact \ +-d api.token=${API_TOKEN} \ +-d buildTargetPHID=${TARGET_PHID} \ +-d artifactKey="Jenkins URI" \ +-d artifactType=uri \ +-d artifactData[uri]=${BUILD_URL} \ +-d artifactData[name]="View Jenkins result" \ +-d artifactData[ui.external]=1 +""" +} diff --git a/cmake/AkantuExtraCompilationProfiles.cmake b/cmake/AkantuExtraCompilationProfiles.cmake index b2e5bc731..36214b7eb 100644 --- a/cmake/AkantuExtraCompilationProfiles.cmake +++ b/cmake/AkantuExtraCompilationProfiles.cmake @@ -1,34 +1,34 @@ #Profiling -set(CMAKE_CXX_FLAGS_PROFILING "-g -pg -DNDEBUG -DAKANTU_NDEBUG -O2" +set(CMAKE_CXX_FLAGS_PROFILING "-g -gdb3 -pg -DNDEBUG -DAKANTU_NDEBUG -O2" CACHE STRING "Flags used by the compiler during profiling builds") -set(CMAKE_C_FLAGS_PROFILING "-g -pg -DNDEBUG -DAKANTU_NDEBUG -O2" +set(CMAKE_C_FLAGS_PROFILING "-g -gdb3 -pg -DNDEBUG -DAKANTU_NDEBUG -O2" CACHE STRING "Flags used by the compiler during profiling builds") -set(CMAKE_Fortran_FLAGS_PROFILING "-g -pg -DNDEBUG -DAKANTU_NDEBUG -O2" +set(CMAKE_Fortran_FLAGS_PROFILING "-g -gdb3 -pg -DNDEBUG -DAKANTU_NDEBUG -O2" CACHE STRING "Flags used by the compiler during profiling builds") set(CMAKE_EXE_LINKER_FLAGS_PROFILING "-pg" CACHE STRING "Flags used by the linker during profiling builds") set(CMAKE_SHARED_LINKER_FLAGS_PROFILING "-pg" CACHE STRING "Flags used by the linker during profiling builds") # Sanitize the code if ((CMAKE_CXX_COMPILER_ID STREQUAL "GNU" AND CMAKE_CXX_COMPILER_VERSION VERSION_GREATER "5.2") OR CMAKE_CXX_COMPILER_ID STREQUAL "Clang") - set(CMAKE_CXX_FLAGS_SANITIZE "-g -O2 -fsanitize=address -fsanitize=leak -fsanitize=undefined" + set(CMAKE_CXX_FLAGS_SANITIZE "-g -gdb3 -O2 -fsanitize=address -fsanitize=leak -fsanitize=undefined -fno-omit-frame-pointer" CACHE STRING "Flags used by the compiler during sanitining builds") - set(CMAKE_C_FLAGS_SANITIZE "-g -O2 -fsanitize=address -fsanitize=leak -fsanitize=undefined" + set(CMAKE_C_FLAGS_SANITIZE "-g -gdb3 -O2 -fsanitize=address -fsanitize=leak -fsanitize=undefined -fno-omit-frame-pointer" CACHE STRING "Flags used by the compiler during sanitining builds") - set(CMAKE_Fortran_FLAGS_SANITIZE "-g -O2 -fsanitize=address -fsanitize=leak -fsanitize=undefined" + set(CMAKE_Fortran_FLAGS_SANITIZE "-g -gdb3 -O2 -fsanitize=address -fsanitize=leak -fsanitize=undefined -fno-omit-frame-pointer" CACHE STRING "Flags used by the compiler during sanitining builds") - set(CMAKE_EXE_LINKER_FLAGS_SANITIZE "-fsanitize=address -fsanitize=leak -fsanitize=undefined" + set(CMAKE_EXE_LINKER_FLAGS_SANITIZE "-fsanitize=address -fsanitize=leak -fsanitize=undefined -fno-omit-frame-pointer" CACHE STRING "Flags used by the linker during sanitining builds") - set(CMAKE_SHARED_LINKER_FLAGS_SANITIZE "-fsanitize=address -fsanitize=leak -fsanitize=undefined" + set(CMAKE_SHARED_LINKER_FLAGS_SANITIZE "-fsanitize=address -fsanitize=leak -fsanitize=undefined -fno-omit-frame-pointer" CACHE STRING "Flags used by the linker during sanitining builds") mark_as_advanced(CMAKE_CXX_FLAGS_PROFILING CMAKE_C_FLAGS_PROFILING CMAKE_Fortran_FLAGS_PROFILING CMAKE_EXE_LINKER_FLAGS_PROFILING CMAKE_SHARED_LINKER_FLAGS_PROFILING CMAKE_SHARED_LINKER_FLAGS_SANITIZE CMAKE_CXX_FLAGS_SANITIZE CMAKE_C_FLAGS_SANITIZE CMAKE_Fortran_FLAGS_SANITIZE CMAKE_EXE_LINKER_FLAGS_SANITIZE ) endif() diff --git a/cmake/akantu_environement.csh.in b/cmake/akantu_environement.csh.in index 55e6493f4..3bd1c5304 100644 --- a/cmake/akantu_environement.csh.in +++ b/cmake/akantu_environement.csh.in @@ -1,2 +1,2 @@ -setenv PYTHONPATH $PYTHONPATH:@PROJECT_BINARY_DIR@/python -setenv LD_LIBRARY_PATH ${LD_LIBRARY_PATH}:@PROJECT_BINARY_DIR@/src +setenv PYTHONPATH $PYTHONPATH:"@PROJECT_BINARY_DIR@/python" +setenv LD_LIBRARY_PATH ${LD_LIBRARY_PATH}:"@PROJECT_BINARY_DIR@/src" diff --git a/cmake/akantu_environement.sh.in b/cmake/akantu_environement.sh.in index b8640c9b3..762163e14 100644 --- a/cmake/akantu_environement.sh.in +++ b/cmake/akantu_environement.sh.in @@ -1,2 +1,2 @@ -export PYTHONPATH=$PYTHONPATH:@PROJECT_BINARY_DIR@/python:@PROJECT_SOURCE_DIR@/test:@PROJECT_SOURCE_DIR@/test/test_fe_engine -export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:@PROJECT_BINARY_DIR@/src +export PYTHONPATH=$PYTHONPATH:"@PROJECT_BINARY_DIR@/python":"@PROJECT_SOURCE_DIR@/test":"@PROJECT_SOURCE_DIR@/test/test_fe_engine" +export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:"@PROJECT_BINARY_DIR@/src" diff --git a/doc/manual/manual-gettingstarted.tex b/doc/manual/manual-gettingstarted.tex index c2912282a..383b48e38 100644 --- a/doc/manual/manual-gettingstarted.tex +++ b/doc/manual/manual-gettingstarted.tex @@ -1,462 +1,466 @@ \chapter{Getting Started} \section{Downloading the Code} The \akantu source code can be requested using the form accessible at the URL \url{http://lsms.epfl.ch/akantu}. There, you will be asked to accept the LGPL license terms. \section{Compiling \akantu} \akantu is a \code{cmake} project, so to configure it, you can either follow the usual way: \begin{command} > cd akantu > mkdir build > cd build > ccmake .. [ Set the options that you need ] > make > make install \end{command} \noindent Or, use the \code{Makefile} we added for your convenience to handle the \code{cmake} configuration \begin{command} > cd akantu > make config > make > make install \end{command} \noindent All the \akantu options are documented in Appendix \ref{app:package-dependencies}. \section{Writing a \texttt{main} Function\label{sect:common:main}} \label{sec:writing_main} First of all, \akantu needs to be initialized. The memory management included in the core library handles the correct allocation and de-allocation of vectors, structures and/or objects. Moreover, in parallel computations, the initialization procedure performs the communication setup. This is achieved by a pair of functions (\code{initialize} and \code{finalize}) that are used as follows: \begin{cpp} #include "aka_common.hh" #include "..." using namespace akantu; int main(int argc, char *argv[]) { initialize("input_file.dat", argc, argv); // your code ... finalize(); } \end{cpp} The \code{initialize} function takes the text inpute file and the program parameters which can be parsed by \akantu in due form (see \ref{sect:parser}). Obviously it is necessary to include all files needed in main. In this manual all provided code implies the usage of \code{akantu} as namespace. \section{Creating and Loading a Mesh\label{sect:common:mesh}} In its current state, \akantu supports three types of meshes: Gmsh~\cite{gmsh}, Abaqus~\cite{abaqus} and Diana~\cite{diana}. Once a \code{Mesh} object is created with a given spatial dimension, it can be filled by reading a mesh input file. The method \code{read} of the class \code{Mesh} infers the mesh type from the file extension. If a non-standard file extension is used, the mesh type has to be specified. \begin{cpp} UInt spatial_dimension = 2; Mesh mesh(spatial_dimension); // Reading Gmsh files mesh.read("my_gmsh_mesh.msh"); mesh.read("my_gmsh_mesh", _miot_gmsh); // Reading Abaqus files mesh.read("my_abaqus_mesh.inp"); mesh.read("my_abaqus_mesh", _miot_abaqus); // Reading Diana files mesh.read("my_diana_mesh.dat"); mesh.read("my_diana_mesh", _miot_diana); \end{cpp} The Gmsh reader adds the geometrical and physical tags as mesh data. The physical values are stored as a \code{UInt} data called \code{tag\_0}, if a string name is provided it is stored as a \code{std::string} data named \code{physical\_names}. The geometrical tag is stored as a \code{UInt} data named \code{tag\_1}. The Abaqus reader stores the \code{ELSET} in ElementGroups and the \code{NSET} in NodeGroups. The material assignment can be retrieved from the \code{std::string} mesh data named \code{abaqus\_material}. % \akantu supports meshes generated with Gmsh~\cite{gmsh}, a free % software available at \url{http://geuz.org/gmsh/} where a detailed % documentation can be found. Consequently, this manual will not provide % Gmsh usage directions. Gmsh outputs meshes in \code{.msh} format that can be read % by \akantu. In order to import a mesh, it is necessary to create % a \code{Mesh} object through the following function calls: % \begin{cpp} % UInt spatial_dimension = 2; % Mesh mesh(spatial_dimension); % \end{cpp} % The only parameter that has to be specified by the user is the spatial % dimension of the problem. Now it is possible to read a \code{.msh} file with % a \code{MeshIOMSH} object that takes care of loading a mesh to memory. % This step is carried out by: % \begin{cpp} % mesh.read("square.msh"); % \end{cpp} % where the \code{MeshIOMSH} object is first created before being % used to read the \code{.msh} file. The mesh file name as well as the \code{Mesh} % object must be specified by the user. % The \code{MeshIOMSH} object can also write mesh files. This feature % is useful to save a mesh that has been modified during a % simulation. The \code{write} method takes care of it: % \begin{cpp} % mesh_io.write("square_modified.msh", mesh); % \end{cpp} % which works exactly like the \code{read} method. % \akantu supports also meshes generated by % DIANA (\url{http://tnodiana.com}), but only in reading mode. A similar % procedure applies where the only % difference is that the \code{MeshIODiana} object should be used % instead of the \code{MeshIOMSH} one. Additional mesh readers can be % introduced into \akantu by coding new \code{MeshIO} classes. \section{Using \texttt{Arrays}} Data in \akantu can be stored in data containers implemented by the \code{Array} class. In its most basic usage, the \code{Array} class implemented in \akantu is similar to the \code{vector} class of the Standard Template Library (STL) for C++. A simple \code{Array} containing a sequence of \code{nb\_element} values (of a given type) can be generated with: \begin{cpp} Array example_array(nb_element); \end{cpp} where \code{type} usually is \code{Real}, \code{Int}, \code{UInt} or \code{bool}. Each value is associated to an index, so that data can be accessed by typing: \begin{cpp} auto & val = example_array(index) \end{cpp} \code{Arrays} can also contain tuples of values for each index. In that case, the number of components per tuple must be specified at the \code{Array} creation. For example, if we want to create an \code{Array} to store the coordinates (sequences of three values) of ten nodes, the appropriate code is the following: \begin{cpp} UInt nb_nodes = 10; UInt spatial_dimension = 3; Array position(nb_nodes, spatial_dimension); \end{cpp} In this case the $x$ position of the eighth node number will be given by \code{position(7, 0)} (in C++, numbering starts at 0 and not 1). If the number of components for the sequences is not specified, the default value of 1 is used. Here is a list of some basic operations that can be performed on \code{Array}: \begin{itemize} \item \code{resize(size)} change the size of the \code{Array}. \item \code{clear()} set all entries of the \code{Array} to zero. \item \code{set(t)} set all entries of the \code{Array} to \code{t}. \item \code{copy(const Array $\&$ other)} copy another \code{Array} into the current one. The two \code{Array} should have the same number of components. \item \code{push$\_$back(tuple)} append a tuple with the correct number of components at the end of the \code{Array}. \item \code{erase(i) erase the value at the i-th position.} \item \code{find(value)} search \code{value} in the current \code{Array}. Return position index of the first occurence or $-1$ if not found. \item \code{storage()} Return the address of the allocated memory of the \code{Array}. \end{itemize} \subsection{\texttt{Arrays} iterators} It is very common in \akantu to loop over arrays to perform a specific treatment. This ranges from geometric calculation on nodal quantities to tensor algebra (in constitutive laws for example). The \code{Array} object has the possibility to request iterators in order to make the writing of loops easier and enhance readability. For instance, a loop over the nodal coordinates can be performed like: \begin{cpp} //accessing the nodal coordinates Array (spatial_dimension components) const auto & nodes = mesh.getNodes(); //creating the iterators auto it = nodes.begin(spatial_dimension); auto end = nodes.end(spatial_dimension); for (; it != end; ++it){ const auto & coords = (*it); //do what you need .... } \end{cpp} In that example, each \code{coords} is a \code{Vector} containing geometrical array of size \code{spatial\_dimension} and the iteration is conveniently performed by the \code{Array} iterator. With the switch to \code{c++14} this can be also written as \begin{cpp} //accessing the nodal coordinates Array (spatial_dimension components) const auto & nodes = mesh.getNodes(); for (const auto & coords : make_view(nodes, spatial_dimension) { //do what you need .... } \end{cpp} The \code{Array} object is intensively used to store second order tensor values. In that case, it should be specified that the returned object type is a matrix when constructing the iterator. This is done when calling the \code{begin} function. For instance, assuming that we have a \code{Array} storing stresses, we can loop over the stored tensors by: \begin{cpp} //creating the iterators auto it = stresses.begin(spatial_dimension, spatial_dimension); auto end = stresses.end(spatial_dimension, spatial_dimension); for (; it != end; ++it){ Matrix & stress = (*it); //do what you need .... } \end{cpp} In that last example, the \code{Matrix} objects are \code{spatial\_dimension} $\times$ \code{spatial\_dimension} matrices. The light objects \code{Matrix} and \code{Vector} can be used and combined to do most common linear algebra. If the number of component is 1, it is possible to use a scalar\_iterator rather than the vector/matrix one. In general, a mesh consists of several kinds of elements. Consequently, the amount of data to be stored can differ for each element type. The straightforward example is the connectivity array, namely the sequences of nodes belonging to each element (linear triangular elements have fewer nodes than, say, rectangular quadratic elements etc.). A particular data structure called \code{ElementTypeMapArray} is provided to easily manage this kind of data. It consists of a group of \code{Arrays}, each associated to an element type. The following code can retrieve the \code{ElementTypeMapArray} which stores the connectivity arrays for a mesh: \begin{cpp} const ElementTypeMapArray & connectivities = mesh.getConnectivities(); \end{cpp} Then, the specific array associated to a given element type can be obtained by \begin{cpp} const Array & connectivity_triangle = connectivities(_triangle_3); \end{cpp} where the first order 3-node triangular element was used in the presented piece of code. \subsection{Vector \& Matrix} The \code{Array} iterators as presented in the previous section can be shaped as \code{Vector} or \code{Matrix}. This objects represent $1^{st}$ and $2^{nd}$ order tensors. As such they come with some functionalities that we will present a bit more into detail in this here. \subsubsection{\texttt{Vector}} \begin{enumerate} \item Accessors: \begin{itemize} \item \code{v(i)} gives the $i^{th}$ component of the vector \code{v} \item \code{v[i]} gives the $i^{th}$ component of the vector \code{v} \item \code{v.size()} gives the number of component \end{itemize} \item Level 1: (results are scalars) \begin{itemize} \item \code{v.norm()} returns the geometrical norm ($L_2$) \item \code{v.norm()} returns the $L_N$ norm defined as $\left(\sum_i |\code{v(i)}|^N\right)^{1/N}$. N can take any positive integer value. There are also some particular values for the most commonly used norms, \code{L\_1} for the Manhattan norm, \code{L\_2} for the geometrical norm and \code{L\_inf} for the norm infinity. \item \code{v.dot(x)} return the dot product of \code{v} and \code{x} \item \code{v.distance(x)} return the geometrical norm of $\code{v} - \code{x}$ \end{itemize} \item Level 2: (results are vectors) \begin{itemize} \item \code{v += s}, \code{v -= s}, \code{v *= s}, \code{v /= s} those are element-wise operators that sum, substract, multiply or divide all the component of \code{v} by the scalar \code{s} \item \code{v += x}, \code{v -= x} sums or substracts the vector \code{x} to/from \code{v} \item \code{v.mul(A, x, alpha)} stores the result of $\alpha \mat{A} \vec{x}$ in \code{v}, $\alpha$ is equal to 1 by default \item \code{v.solve(A, b)} stores the result of the resolution of the system $\mat{A} \vec{x} = \vec{b}$ in \code{v} \item \code{v.crossProduct(v1, v2)} computes the cross product of \code{v1} and \code{v2} and stores the result in \code{v} \end{itemize} \end{enumerate} \subsubsection{\texttt{Matrix}} \begin{enumerate} \item Accessors: \begin{itemize} \item \code{A(i, j)} gives the component $A_{ij}$ of the matrix \code{A} \item \code{A(i)} gives the $i^{th}$ column of the matrix as a \code{Vector} \item \code{A[k]} gives the $k^{th}$ component of the matrix, matrices are stored in a column major way, which means that to access $A_{ij}$, $k = i + j M$ \item \code{A.rows()} gives the number of rows of \code{A} ($M$) \item \code{A.cols()} gives the number of columns of \code{A} ($N$) \item \code{A.size()} gives the number of component in the matrix ($M \times N$) \end{itemize} \item Level 1: (results are scalars) \begin{itemize} \item \code{A.norm()} is equivalent to \code{A.norm()} \item \code{A.norm()} returns the $L_N$ norm defined as $\left(\sum_i\sum_j |\code{A(i,j)}|^N\right)^{1/N}$. N can take any positive integer value. There are also some particular values for the most commonly used norms, \code{L\_1} for the Manhattan norm, \code{L\_2} for the geometrical norm and \code{L\_inf} for the norm infinity. \item \code{A.trace()} return the trace of \code{A} \item \code{A.det()} return the determinant of \code{A} \item \code{A.doubleDot(B)} return the double dot product of \code{A} and \code{B}, $\mat{A}:\mat{B}$ \end{itemize} \item Level 3: (results are matrices) \begin{itemize} \item \code{A.eye(s)}, \code{Matrix::eye(s)} fills/creates a matrix with the $s\mat{I}$ with $\mat{I}$ the identity matrix \item \code{A.inverse(B)} stores $\mat{B}^{-1}$ in \code{A} \item \code{A.transpose()} returns $\mat{A}^{t}$ \item \code{A.outerProduct(v1, v2)} stores $\vec{v_1} \vec{v_2}^{t}$ in \code{A} \item \code{C.mul(A, B, alpha)}: stores the result of the product of \code{A} and code{B} time the scalar \code{alpha} in \code{C}. \code{t\_A} and \code{t\_B} are boolean defining if \code{A} and \code{B} should be transposed or not. \begin{tabular}{ccl} \toprule \code{t\_A} & \code{t\_B} & result \\ \midrule false & false & $\mat{C} = \alpha \mat{A} \mat{B}$\\ false & true & $\mat{C} = \alpha \mat{A} \mat{B}^t$\\ true & false & $\mat{C} = \alpha \mat{A}^t \mat{B}$\\ true & true & $\mat{C} = \alpha \mat{A}^t \mat{B}^t$\\ \bottomrule \end{tabular} \item \code{A.eigs(d, V)} this method computes the eigenvalues and eigenvectors of \code{A} and store the results in \code{d} and \code{V} such that $\code{d(i)} = \lambda_i$ and $\code{V(i)} = \vec{v_i}$ with $\mat{A}\vec{v_i} = \lambda_i\vec{v_i}$ and $\lambda_1 > ... > \lambda_i > ... > \lambda_N$ \end{itemize} \end{enumerate} \subsubsection{\texttt{Tensor3}} Accessors: \begin{itemize} \item \code{t(i, j, k)} gives the component $T_{ijk}$ of the tensor \code{t} \item \code{t(k)} gives the $k^{th}$ two-dimensional tensor as a \code{Matrix} \item \code{t[k]} gives the $k^{th}$ two-dimensional tensor as a \code{Matrix} \end{itemize} \section{Manipulating group of nodes and/or elements\label{sect:common:groups}} \akantu provides the possibility to manipulate subgroups of elements and nodes. Any \code{ElementGroup} and/or \code{NodeGroup} must be managed by a \code{GroupManager}. Such a manager has the role to associate group objects to names. This is a useful feature, in particular for the application of the boundary conditions, as will be demonstrated in section \ref{sect:smm:boundary}. To most general group manager is the \code{Mesh} class which inheritates from the \code{GroupManager} class. For instance, the following code shows how to request an element group to a mesh: \begin{cpp} // request creation of a group of nodes NodeGroup & my_node_group = mesh.createNodeGroup("my_node_group"); // request creation of a group of elements ElementGroup & my_element_group = mesh.createElementGroup("my_element_group"); /* fill and use the groups */ \end{cpp} \subsection{The \texttt{NodeGroup} object} A group of nodes is stored in \code{NodeGroup} objects. They are quite simple objects which store the indexes of the selected nodes in a \code{Array}. Nodes are selected by adding them when calling \code{NodeGroup::add}. For instance you can select nodes having a positive $X$ coordinate with the following code: \begin{cpp} const auto & nodes = mesh.getNodes(); auto & group = mesh.createNodeGroup("XpositiveNode"); for (auto && data : enumerate(make_view(nodes, spatial_dimension))){ auto node = std::get<0>(data); const auto & position = std::get<1>(data); if (position(0) > 0) group.add(node); } \end{cpp} \subsection{The \texttt{ElementGroup} object} A group of elements is stored in \code{ElementGroup} objects. Since a group can contain elements of various types the \code{ElementGroup} object stores indexes in a \code{ElementTypeMapArray} object. Then elements can be added to the group by calling \code{addElement}. For instance, selecting the elements for which the barycenter of the nodes has a positive $X$ coordinate can be made with: \begin{cpp} auto & group = mesh.createElementGroup("XpositiveElement"); Vector barycenter(spatial_dimension); for(auto type : mesh.elementTypes()){ UInt nb_element = mesh.getNbElement(type); for(UInt e = 0; e < nb_element; ++e) { Element element{type, e, _not_ghost}; mesh.getBarycenter(element, barycenter); if (barycenter(_x) > 0.) group.add(element); } } \end{cpp} \section{Compiling your simulation} -The easiest way it to create a \code{cmake} project. The minimum \code{CMakeLists.txt} file would look like +The easiest way to compile your simulation is to create a \code{cmake} project by putting all your code in some directory of your choosing. Then, make sure that you have \code{cmake} installed and create a \code{CMakeLists.txt} file. An example of a minimal \code{CMakeLists.txt} file would look like this: \begin{cmake} project(my_simu) cmake_minimum_required(VERSION 3.0.0) find_package(Akantu REQUIRED) - + add_akantu_simulation(my_simu my_simu.cc) \end{cmake} +% +Then create a directory called \code{build} and inside it execute \code{ccmake ..} which opens a configuration screen. If you installed \akantu in a standard directory such as \code{/usr/local} (using \code{make install}), you should be able to compile by hitting the \code{c} key, setting \code{CMAKE\textunderscore{}BUILD\textunderscore{}TYPE} to \code{Release}, hitting the \code{c} key again a few times and then finishing the configuration with the \code{g} key. After that, you can type \code{make} to build your simulation. If you change your simulation code later on, you only need to type \code{make} again. + +If you get an error that \code{FindAkantu.cmake} was not found, you have to set the \code{Akantu\textunderscore{}DIR} variable, which will appear after dismissing the error message. If you built \akantu without running \code{make install}, the variable should be set to the \code{build} subdirectory of the \akantu source code. If you installed it in \code{\$PREFIX}, set the variable to \code{\$PREFIX/share/cmake/Akantu} instead. %%% Local Variables: %%% mode: latex %%% TeX-master: "manual" %%% End: diff --git a/examples/explicit/explicit_dynamic.cc b/examples/explicit/explicit_dynamic.cc index e640a7344..c03506e49 100644 --- a/examples/explicit/explicit_dynamic.cc +++ b/examples/explicit/explicit_dynamic.cc @@ -1,106 +1,113 @@ /** * @file explicit_dynamic.cc * * @author Seyedeh Mohadeseh Taheri Mousavi * * @date creation: Sun Jul 12 2015 * @date last modification: Mon Jan 18 2016 * * @brief This code refers to the explicit dynamic example from the user manual * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ using namespace akantu; int main(int argc, char * argv[]) { initialize("material.dat", argc, argv); const UInt spatial_dimension = 3; const Real pulse_width = 2.; const Real A = 0.01; Real time_step; Real time_factor = 0.8; UInt max_steps = 1000; Mesh mesh(spatial_dimension); - mesh.read("bar.msh"); + if(Communicator::getStaticCommunicator().whoAmI() == 0) + mesh.read("bar.msh"); + + mesh.distribute(); + + mesh.makePeriodic(_x); + mesh.makePeriodic(_y); + mesh.makePeriodic(_z); SolidMechanicsModel model(mesh); /// model initialization model.initFull(_analysis_method = _explicit_lumped_mass); time_step = model.getStableTimeStep(); std::cout << "Time Step = " << time_step * time_factor << "s (" << time_step << "s)" << std::endl; model.setTimeStep(time_step * time_factor); /// boundary and initial conditions Array & displacement = model.getDisplacement(); const Array & nodes = mesh.getNodes(); for (UInt n = 0; n < mesh.getNbNodes(); ++n) { - Real x = nodes(n); + Real x = nodes(n) - 2; // Sinus * Gaussian Real L = pulse_width; Real k = 0.1 * 2 * M_PI * 3 / L; displacement(n) = A * sin(k * x) * exp(-(k * x) * (k * x) / (L * L)); } std::ofstream energy; energy.open("energy.csv"); energy << "id,rtime,epot,ekin,tot" << std::endl; model.setBaseName("explicit_dynamic"); model.addDumpField("displacement"); model.addDumpField("velocity"); model.addDumpField("acceleration"); model.addDumpField("stress"); model.dump(); for (UInt s = 1; s <= max_steps; ++s) { model.solveStep(); Real epot = model.getEnergy("potential"); Real ekin = model.getEnergy("kinetic"); energy << s << "," << s * time_step << "," << epot << "," << ekin << "," << epot + ekin << "," << std::endl; if (s % 10 == 0) std::cout << "passing step " << s << "/" << max_steps << std::endl; model.dump(); } energy.close(); finalize(); return EXIT_SUCCESS; } diff --git a/extra_packages/igfem/src/integrator_gauss_igfem_inline_impl.cc b/extra_packages/igfem/src/integrator_gauss_igfem_inline_impl.cc index c027b3a9a..8f75d8395 100644 --- a/extra_packages/igfem/src/integrator_gauss_igfem_inline_impl.cc +++ b/extra_packages/igfem/src/integrator_gauss_igfem_inline_impl.cc @@ -1,452 +1,452 @@ /** * @file integrator_gauss_igfem.hh * * @author Aurelia Isabel Cuba Ramos * * * @brief Inline functions of gauss integrator for the case of IGFEM * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ __END_AKANTU__ #include "fe_engine.hh" #if defined(AKANTU_DEBUG_TOOLS) #include "aka_debug_tools.hh" #endif __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ #define INIT_INTEGRATOR(type) \ computeQuadraturePoints(ghost_type); \ precomputeJacobiansOnQuadraturePoints(nodes, ghost_type); \ checkJacobians(ghost_type); template inline void IntegratorGauss<_ek_igfem, IOF>::initIntegrator(const Array & nodes, const ElementType & type, const GhostType & ghost_type) { AKANTU_BOOST_IGFEM_ELEMENT_SWITCH(INIT_INTEGRATOR); } #undef INIT_INTEGRATOR /* -------------------------------------------------------------------------- */ template template inline UInt IntegratorGauss<_ek_igfem, IOF>::getNbIntegrationPoints( const GhostType &) const { const ElementType sub_type_1 = ElementClassProperty::sub_element_type_1; const ElementType sub_type_2 = ElementClassProperty::sub_element_type_2; UInt nb_quad_points_sub_1 = GaussIntegrationElement::getNbQuadraturePoints(); UInt nb_quad_points_sub_2 = GaussIntegrationElement::getNbQuadraturePoints(); return (nb_quad_points_sub_1 + nb_quad_points_sub_2); } /* -------------------------------------------------------------------------- */ template template inline void IntegratorGauss<_ek_igfem, IOF>::integrateOnElement( const Array & f, Real * intf, UInt nb_degree_of_freedom, const UInt elem, const GhostType & ghost_type) const { Array & jac_loc = jacobians(type, ghost_type); UInt nb_quadrature_points = getNbIntegrationPoints(); AKANTU_DEBUG_ASSERT(f.getNbComponent() == nb_degree_of_freedom, "The vector f do not have the good number of component."); Real * f_val = f.storage() + elem * f.getNbComponent(); Real * jac_val = jac_loc.storage() + elem * nb_quadrature_points; integrate(f_val, jac_val, intf, nb_degree_of_freedom, nb_quadrature_points); } /* -------------------------------------------------------------------------- */ template template inline Real IntegratorGauss<_ek_igfem, IOF>::integrate( const Vector & in_f, UInt index, const GhostType & ghost_type) const { const Array & jac_loc = jacobians(type, ghost_type); UInt nb_quadrature_points = getNbIntegrationPoints(); AKANTU_DEBUG_ASSERT(in_f.size() == nb_quadrature_points, "The vector f do not have nb_quadrature_points entries."); Real * jac_val = jac_loc.storage() + index * nb_quadrature_points; Real intf; integrate(in_f.storage(), jac_val, &intf, 1, nb_quadrature_points); return intf; return 0.; } /* -------------------------------------------------------------------------- */ template inline void IntegratorGauss<_ek_igfem, IOF>::integrate(Real * f, Real * jac, Real * inte, UInt nb_degree_of_freedom, UInt nb_quadrature_points) const { memset(inte, 0, nb_degree_of_freedom * sizeof(Real)); Real * cjac = jac; for (UInt q = 0; q < nb_quadrature_points; ++q) { for (UInt dof = 0; dof < nb_degree_of_freedom; ++dof) { inte[dof] += *f * *cjac; ++f; } ++cjac; } } /* -------------------------------------------------------------------------- */ template template inline const Matrix & IntegratorGauss<_ek_igfem, IOF>::getIntegrationPoints( const GhostType & ghost_type) const { AKANTU_DEBUG_ASSERT( quadrature_points.exists(type, ghost_type), "Quadrature points for type " << quadrature_points.printType(type, ghost_type) << " have not been initialized." << " Did you use 'computeQuadraturePoints' function ?"); return quadrature_points(type, ghost_type); } /* -------------------------------------------------------------------------- */ template template inline void IntegratorGauss<_ek_igfem, IOF>::computeQuadraturePoints( const GhostType & ghost_type) { /// typedef for the two subelement_types and the parent element type const ElementType sub_type_1 = ElementClassProperty::sub_element_type_1; const ElementType sub_type_2 = ElementClassProperty::sub_element_type_2; /// store the quadrature points on the two subelements Matrix & quads_sub_1 = quadrature_points(sub_type_1, ghost_type); Matrix & quads_sub_2 = quadrature_points(sub_type_2, ghost_type); quads_sub_1 = GaussIntegrationElement::getQuadraturePoints(); quads_sub_2 = GaussIntegrationElement::getQuadraturePoints(); /// store all quad points for the current type UInt nb_quad_points_sub_1 = GaussIntegrationElement::getNbQuadraturePoints(); UInt nb_quad_points_sub_2 = GaussIntegrationElement::getNbQuadraturePoints(); UInt spatial_dimension = mesh.getSpatialDimension(); Matrix & quads = quadrature_points(type, ghost_type); quads = Matrix(spatial_dimension, nb_quad_points_sub_1 + nb_quad_points_sub_2); Matrix quads_1(quads.storage(), quads.rows(), nb_quad_points_sub_1); quads_1 = quads_sub_1; Matrix quads_2(quads.storage() + quads.rows() * nb_quad_points_sub_1, quads.rows(), nb_quad_points_sub_2); quads_2 = quads_sub_2; } /* -------------------------------------------------------------------------- */ template template inline void IntegratorGauss<_ek_igfem, IOF>::computeJacobianOnQuadPointsByElement( const Matrix & node_coords, Vector & jacobians) { /// optimize: get the matrix from the ElementTypeMap Matrix quad = GaussIntegrationElement::getQuadraturePoints(); // jacobian ElementClass::computeJacobian(quad, node_coords, jacobians); } /* -------------------------------------------------------------------------- */ template inline IntegratorGauss<_ek_igfem, IOF>::IntegratorGauss( const Mesh & mesh, const ID & id, const MemoryID & memory_id) : Integrator(mesh, id, memory_id) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template inline void IntegratorGauss<_ek_igfem, IOF>::checkJacobians( const GhostType & ghost_type) const { AKANTU_DEBUG_IN(); /// typedef for the two subelement_types and the parent element type const ElementType sub_type_1 = ElementClassProperty::sub_element_type_1; const ElementType sub_type_2 = ElementClassProperty::sub_element_type_2; UInt nb_quad_points_sub_1 = GaussIntegrationElement::getNbQuadraturePoints(); UInt nb_quad_points_sub_2 = GaussIntegrationElement::getNbQuadraturePoints(); UInt nb_quadrature_points = nb_quad_points_sub_1 + nb_quad_points_sub_2; UInt nb_element; nb_element = mesh.getConnectivity(type, ghost_type).getSize(); Real * jacobians_val = jacobians(type, ghost_type).storage(); for (UInt i = 0; i < nb_element * nb_quadrature_points; ++i, ++jacobians_val) { if (*jacobians_val < 0) - AKANTU_DEBUG_ERROR( + AKANTU_ERROR( "Negative jacobian computed," << " possible problem in the element node ordering (Quadrature Point " << i % nb_quadrature_points << ":" << i / nb_quadrature_points << ":" << type << ":" << ghost_type << ")"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template inline void IntegratorGauss<_ek_igfem, IOF>::precomputeJacobiansOnQuadraturePoints( const Array & nodes, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); /// typedef for the two subelement_types and the parent element type const ElementType sub_type_1 = ElementClassProperty::sub_element_type_1; const ElementType sub_type_2 = ElementClassProperty::sub_element_type_2; UInt spatial_dimension = mesh.getSpatialDimension(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); /// get the number of nodes for the subelements and the parent element UInt nb_nodes_sub_1 = ElementClass::getNbNodesPerInterpolationElement(); UInt nb_nodes_sub_2 = ElementClass::getNbNodesPerInterpolationElement(); UInt nb_quadrature_points_sub_1 = GaussIntegrationElement::getNbQuadraturePoints(); UInt nb_quadrature_points_sub_2 = GaussIntegrationElement::getNbQuadraturePoints(); UInt nb_quadrature_points = nb_quadrature_points_sub_1 + nb_quadrature_points_sub_2; UInt nb_element = mesh.getNbElement(type, ghost_type); Array * jacobians_tmp; if (!jacobians.exists(type, ghost_type)) jacobians_tmp = &jacobians.alloc(nb_element * nb_quadrature_points, 1, type, ghost_type); else { jacobians_tmp = &jacobians(type, ghost_type); jacobians_tmp->resize(nb_element * nb_quadrature_points); } Array::vector_iterator jacobians_it = jacobians_tmp->begin_reinterpret(nb_quadrature_points, nb_element); Vector weights_sub_1 = GaussIntegrationElement::getWeights(); Vector weights_sub_2 = GaussIntegrationElement::getWeights(); Array x_el(0, spatial_dimension * nb_nodes_per_element); FEEngine::extractNodalToElementField(mesh, nodes, x_el, type, ghost_type); Array::const_matrix_iterator x_it = x_el.begin(spatial_dimension, nb_nodes_per_element); // Matrix local_coord(spatial_dimension, nb_nodes_per_element); for (UInt elem = 0; elem < nb_element; ++elem, ++jacobians_it, ++x_it) { const Matrix & X = *x_it; Matrix sub_1_coords(spatial_dimension, nb_nodes_sub_1); Matrix sub_2_coords(spatial_dimension, nb_nodes_sub_2); ElementClass::getSubElementCoords(X, sub_1_coords, 0); ElementClass::getSubElementCoords(X, sub_2_coords, 1); Vector & J = *jacobians_it; /// initialize vectors to store the jacobians for each subelement Vector J_sub_1(nb_quadrature_points_sub_1); Vector J_sub_2(nb_quadrature_points_sub_2); computeJacobianOnQuadPointsByElement(sub_1_coords, J_sub_1); computeJacobianOnQuadPointsByElement(sub_2_coords, J_sub_2); J_sub_1 *= weights_sub_1; J_sub_2 *= weights_sub_2; /// copy results into the jacobian vector for this element for (UInt i = 0; i < nb_quadrature_points_sub_1; ++i) { J(i) = J_sub_1(i); } for (UInt i = 0; i < nb_quadrature_points_sub_2; ++i) { J(i + nb_quadrature_points_sub_1) = J_sub_2(i); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template inline void IntegratorGauss<_ek_igfem, IOF>::integrate( const Array & in_f, Array & intf, UInt nb_degree_of_freedom, const GhostType & ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(jacobians.exists(type, ghost_type), "No jacobians for the type " << jacobians.printType(type, ghost_type)); const Matrix & quads = quadrature_points(type, ghost_type); UInt nb_points = quads.cols(); const Array & jac_loc = jacobians(type, ghost_type); Array::const_matrix_iterator J_it; Array::matrix_iterator inte_it; Array::const_matrix_iterator f_it; UInt nb_element; Array * filtered_J = NULL; if (filter_elements != empty_filter) { nb_element = filter_elements.getSize(); filtered_J = new Array(0, jac_loc.getNbComponent()); FEEngine::filterElementalData(mesh, jac_loc, *filtered_J, type, ghost_type, filter_elements); const Array & cfiltered_J = *filtered_J; // \todo temporary patch J_it = cfiltered_J.begin_reinterpret(nb_points, 1, nb_element); } else { nb_element = mesh.getNbElement(type, ghost_type); J_it = jac_loc.begin_reinterpret(nb_points, 1, nb_element); } intf.resize(nb_element); f_it = in_f.begin_reinterpret(nb_degree_of_freedom, nb_points, nb_element); inte_it = intf.begin_reinterpret(nb_degree_of_freedom, 1, nb_element); for (UInt el = 0; el < nb_element; ++el, ++J_it, ++f_it, ++inte_it) { const Matrix & f = *f_it; const Matrix & J = *J_it; Matrix & inte_f = *inte_it; inte_f.mul(f, J); } delete filtered_J; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template template inline Real IntegratorGauss<_ek_igfem, IOF>::integrate( const Array & in_f, const GhostType & ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(jacobians.exists(type, ghost_type), "No jacobians for the type " << jacobians.printType(type, ghost_type)); Array intfv(0, 1); integrate(in_f, intfv, 1, ghost_type, filter_elements); UInt nb_values = intfv.getSize(); if (nb_values == 0) return 0.; UInt nb_values_to_sum = nb_values >> 1; std::sort(intfv.begin(), intfv.end()); // as long as the half is not empty while (nb_values_to_sum) { UInt remaining = (nb_values - 2 * nb_values_to_sum); if (remaining) intfv(nb_values - 2) += intfv(nb_values - 1); // sum to consecutive values and store the sum in the first half for (UInt i = 0; i < nb_values_to_sum; ++i) { intfv(i) = intfv(2 * i) + intfv(2 * i + 1); } nb_values = nb_values_to_sum; nb_values_to_sum >>= 1; } AKANTU_DEBUG_OUT(); return intfv(0); } /* -------------------------------------------------------------------------- */ template template inline void IntegratorGauss<_ek_igfem, IOF>::integrateOnIntegrationPoints( const Array & in_f, Array & intf, UInt nb_degree_of_freedom, const GhostType & ghost_type, const Array & filter_elements) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(jacobians.exists(type, ghost_type), "No jacobians for the type " << jacobians.printType(type, ghost_type)); UInt nb_element; const Matrix & quads = quadrature_points(type, ghost_type); UInt nb_points = quads.cols(); const Array & jac_loc = jacobians(type, ghost_type); Array::const_scalar_iterator J_it; Array::vector_iterator inte_it; Array::const_vector_iterator f_it; Array * filtered_J = NULL; if (filter_elements != empty_filter) { nb_element = filter_elements.getSize(); filtered_J = new Array(0, jac_loc.getNbComponent()); FEEngine::filterElementalData(mesh, jac_loc, *filtered_J, type, ghost_type, filter_elements); J_it = filtered_J->begin(); } else { nb_element = mesh.getNbElement(type, ghost_type); J_it = jac_loc.begin(); } intf.resize(nb_element * nb_points); f_it = in_f.begin(nb_degree_of_freedom); inte_it = intf.begin(nb_degree_of_freedom); for (UInt el = 0; el < nb_element; ++el, ++J_it, ++f_it, ++inte_it) { const Real & J = *J_it; const Vector & f = *f_it; Vector & inte_f = *inte_it; inte_f = f; inte_f *= J; } delete filtered_J; AKANTU_DEBUG_OUT(); } diff --git a/extra_packages/igfem/src/material_igfem/material_igfem_elastic.cc b/extra_packages/igfem/src/material_igfem/material_igfem_elastic.cc index 20f8ef6d2..e6f6f2c3c 100644 --- a/extra_packages/igfem/src/material_igfem/material_igfem_elastic.cc +++ b/extra_packages/igfem/src/material_igfem/material_igfem_elastic.cc @@ -1,243 +1,243 @@ /** * @file material_igfem_elastic.cc * * @author Aurelia Isabel Cuba Ramos * * * @brief Specializaton of material class for the igfem elastic material * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * */ /* -------------------------------------------------------------------------- */ #include "material_igfem_elastic.hh" #include "material_elastic.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialIGFEMElastic::MaterialIGFEMElastic(SolidMechanicsModel & model, const ID & id) : Material(model, id), Parent(model, id), lambda("lambda", *this), mu("mu", *this), kpa("kappa", *this) { AKANTU_DEBUG_IN(); this->initialize(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialIGFEMElastic::initialize() { this->lambda.initialize(1); this->mu.initialize(1); this->kpa.initialize(1); } /* -------------------------------------------------------------------------- */ template void MaterialIGFEMElastic::initMaterial() { AKANTU_DEBUG_IN(); Parent::initMaterial(); /// insert the sub_material names into the map this->sub_material_names[0] = this->name_sub_mat_1; this->sub_material_names[1] = this->name_sub_mat_2; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialIGFEMElastic::updateElasticInternals( const Array & element_list) { /// compute the Lamé constants for both sub-materials Vector lambda_per_sub_mat(this->nb_sub_materials); Vector mu_per_sub_mat(this->nb_sub_materials); Vector kpa_per_sub_mat(this->nb_sub_materials); for (UInt i = 0; i < this->nb_sub_materials; ++i) { ID mat_name = this->sub_material_names[i]; const MaterialElastic & mat = dynamic_cast &>( this->model->getMaterial(mat_name)); lambda_per_sub_mat(i) = mat.getLambda(); mu_per_sub_mat(i) = mat.getMu(); kpa_per_sub_mat(i) = mat.getKappa(); } for (ghost_type_t::iterator g = ghost_type_t::begin(); g != ghost_type_t::end(); ++g) { GhostType ghost_type = *g; /// loop over all types in the material typedef ElementTypeMapArray::type_iterator iterator; iterator it = this->element_filter.firstType(spatial_dimension, ghost_type, _ek_igfem); iterator last_type = this->element_filter.lastType(spatial_dimension, ghost_type, _ek_igfem); /// loop over all types in the filter for (; it != last_type; ++it) { ElementType el_type = *it; if (el_type == _igfem_triangle_4) this->template setSubMaterial<_igfem_triangle_4>(element_list, ghost_type); else if (el_type == _igfem_triangle_5) this->template setSubMaterial<_igfem_triangle_5>(element_list, ghost_type); else - AKANTU_DEBUG_ERROR( + AKANTU_ERROR( "There is currently no other IGFEM type implemented"); UInt nb_element = this->element_filter(el_type, ghost_type).getSize(); UInt nb_quads = this->fem->getNbIntegrationPoints(el_type); /// get pointer to internals for given type Real * lambda_ptr = this->lambda(el_type, ghost_type).storage(); Real * mu_ptr = this->mu(el_type, ghost_type).storage(); Real * kpa_ptr = this->kpa(el_type, ghost_type).storage(); UInt * sub_mat_ptr = this->sub_material(el_type, ghost_type).storage(); for (UInt q = 0; q < nb_element * nb_quads; ++q, ++lambda_ptr, ++mu_ptr, ++kpa_ptr, ++sub_mat_ptr) { UInt index = *sub_mat_ptr; *lambda_ptr = lambda_per_sub_mat(index); *mu_ptr = mu_per_sub_mat(index); *kpa_ptr = kpa_per_sub_mat(index); } } } } /* -------------------------------------------------------------------------- */ template void MaterialIGFEMElastic::computeStress( ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Parent::computeStress(el_type, ghost_type); if (!this->finite_deformation) { /// get pointer to internals Real * lambda_ptr = this->lambda(el_type, ghost_type).storage(); Real * mu_ptr = this->mu(el_type, ghost_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); this->computeStressOnQuad(grad_u, sigma, *lambda_ptr, *mu_ptr); ++lambda_ptr; ++mu_ptr; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; } else { AKANTU_DEBUG_TO_IMPLEMENT(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialIGFEMElastic::computeTangentModuli( __attribute__((unused)) const ElementType & el_type, Array & tangent_matrix, __attribute__((unused)) GhostType ghost_type) { AKANTU_DEBUG_IN(); /// get pointer to internals Real * lambda_ptr = this->lambda(el_type, ghost_type).storage(); Real * mu_ptr = this->mu(el_type, ghost_type).storage(); MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix); this->computeTangentModuliOnQuad(tangent, *lambda_ptr, *mu_ptr); ++lambda_ptr; ++mu_ptr; MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialIGFEMElastic::computePotentialEnergy( ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); // MaterialThermal::computePotentialEnergy(el_type, // ghost_type); // if(ghost_type != _not_ghost) return; // Array::scalar_iterator epot = this->potential_energy(el_type, // ghost_type).begin(); // if (!this->finite_deformation) { // MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); // this->computePotentialEnergyOnQuad(grad_u, sigma, *epot); // ++epot; // MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; // } else { // Matrix E(spatial_dimension, spatial_dimension); // MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); // this->template gradUToGreenStrain(grad_u, E); // this->computePotentialEnergyOnQuad(E, sigma, *epot); // ++epot; // MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; // } AKANTU_DEBUG_TO_IMPLEMENT(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialIGFEMElastic::computePotentialEnergyByElement( ElementType type, UInt index, Vector & epot_on_quad_points) { // Array::matrix_iterator gradu_it = // this->gradu(type).begin(spatial_dimension, // spatial_dimension); // Array::matrix_iterator gradu_end = // this->gradu(type).begin(spatial_dimension, // spatial_dimension); // Array::matrix_iterator stress_it = // this->stress(type).begin(spatial_dimension, // spatial_dimension); // if (this->finite_deformation) // stress_it = this->piola_kirchhoff_2(type).begin(spatial_dimension, // spatial_dimension); // UInt nb_quadrature_points = // this->model->getFEEngine().getNbQuadraturePoints(type); // gradu_it += index*nb_quadrature_points; // gradu_end += (index+1)*nb_quadrature_points; // stress_it += index*nb_quadrature_points; // Real * epot_quad = epot_on_quad_points.storage(); // Matrix grad_u(spatial_dimension, spatial_dimension); // for(;gradu_it != gradu_end; ++gradu_it, ++stress_it, ++epot_quad) { // if (this->finite_deformation) // this->template gradUToGreenStrain(*gradu_it, // grad_u); // else // grad_u.copy(*gradu_it); // this->computePotentialEnergyOnQuad(grad_u, *stress_it, *epot_quad); // } AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ template void MaterialIGFEMElastic::onElementsAdded( const Array & element_list, const NewElementsEvent & event) { Parent::onElementsAdded(element_list, event); updateElasticInternals(element_list); }; INSTANTIATE_MATERIAL(MaterialIGFEMElastic); __END_AKANTU__ diff --git a/extra_packages/igfem/src/mesh_igfem_spherical_growing_gel_tmpl.hh b/extra_packages/igfem/src/mesh_igfem_spherical_growing_gel_tmpl.hh index c01636c46..4f68e9ad5 100644 --- a/extra_packages/igfem/src/mesh_igfem_spherical_growing_gel_tmpl.hh +++ b/extra_packages/igfem/src/mesh_igfem_spherical_growing_gel_tmpl.hh @@ -1,532 +1,532 @@ /** * @file mesh_igfem_spherical_growing_gel.cc * @author Marco Vocialta * @date Mon Sep 21 12:42:11 2015 * * @brief Implementation of the functions of MeshIgfemSphericalGrowingGel * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "mesh_utils.hh" #include __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MeshIgfemSphericalGrowingGel::MeshIgfemSphericalGrowingGel(Mesh & mesh) : mesh(mesh), nb_nodes_fem(mesh.getNbNodes()), nb_enriched_nodes(0), synchronizer(NULL) {} /* -------------------------------------------------------------------------- */ template void MeshIgfemSphericalGrowingGel::init() { nb_nodes_fem = mesh.getNbNodes(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; Mesh::type_iterator tit = mesh.firstType(dim, ghost_type); Mesh::type_iterator tend = mesh.lastType(dim, ghost_type); for (; tit != tend; ++tit) { // loop to add corresponding IGFEM element types if (*tit == _triangle_3) { mesh.addConnectivityType(_igfem_triangle_4, ghost_type); mesh.addConnectivityType(_igfem_triangle_5, ghost_type); } else - AKANTU_DEBUG_ERROR("Not ready for mesh type " << *tit); + AKANTU_ERROR("Not ready for mesh type " << *tit); } tit = mesh.firstType(dim, ghost_type, _ek_not_defined); tend = mesh.lastType(dim, ghost_type, _ek_not_defined); for (; tit != tend; ++tit) { AKANTU_BOOST_LIST_SWITCH(INSTANTIATOR, ELEMENT_LIST, *tit); } } } /* -------------------------------------------------------------------------- */ template void MeshIgfemSphericalGrowingGel::computeMeshQueryListIntersectionPoint( const std::list & query_list) { /// store number of currently enriched nodes this->nb_enriched_nodes = mesh.getNbNodes() - nb_nodes_fem; UInt nb_old_nodes = mesh.getNbNodes(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; Mesh::type_iterator iit = mesh.firstType(dim, ghost_type, _ek_not_defined); Mesh::type_iterator iend = mesh.lastType(dim, ghost_type, _ek_not_defined); for (; iit != iend; ++iit) { MeshAbstractIntersector & intersector = *intersectors(*iit, ghost_type); intersector.constructData(ghost_type); intersector.computeMeshQueryListIntersectionPoint(query_list, nb_old_nodes); const Array intersection_points_current_type = *(intersector.getIntersectionPoints()); const Array & new_node_per_elem = intersector.getNewNodePerElem(); /// Send the new node event UInt new_points = intersection_points_current_type.getSize(); StaticCommunicator::getStaticCommunicator().allReduce(&new_points, 1, _so_sum); if (new_points > 0) { Array & nodes = this->mesh.getNodes(); UInt nb_node = nodes.getSize(); Array::const_vector_iterator intersec_p_it = intersection_points_current_type.begin(dim); NewIGFEMNodesEvent new_nodes_event; for (UInt in = 0; in < intersection_points_current_type.getSize(); ++in, ++intersec_p_it) { nodes.push_back(*intersec_p_it); new_nodes_event.getList().push_back(nb_node + in); } new_nodes_event.setNewNodePerElem(new_node_per_elem); new_nodes_event.setType(*iit); new_nodes_event.setGhostType(ghost_type); this->mesh.sendEvent(new_nodes_event); } } } } /* -------------------------------------------------------------------------- */ template void MeshIgfemSphericalGrowingGel::removeAdditionalNodes() { AKANTU_DEBUG_IN(); UInt total_nodes = this->mesh.getNbNodes(); UInt nb_new_enriched_nodes = total_nodes - this->nb_enriched_nodes - this->nb_nodes_fem; UInt old_total_nodes = this->nb_nodes_fem + this->nb_enriched_nodes; UInt total_new_nodes = nb_new_enriched_nodes; StaticCommunicator::getStaticCommunicator().allReduce(&total_new_nodes, 1, _so_sum); if (total_new_nodes == 0) return; RemovedNodesEvent remove_nodes(this->mesh); Array & nodes_removed = remove_nodes.getList(); Array & new_numbering = remove_nodes.getNewNumbering(); for (UInt nnod = 0; nnod < this->nb_nodes_fem; ++nnod) { new_numbering(nnod) = nnod; } for (UInt nnod = nb_nodes_fem; nnod < old_total_nodes; ++nnod) { new_numbering(nnod) = UInt(-1); nodes_removed.push_back(nnod); } for (UInt nnod = 0; nnod < nb_new_enriched_nodes; ++nnod) { new_numbering(nnod + old_total_nodes) = this->nb_nodes_fem + nnod; } if (nb_new_enriched_nodes > 0) { for (UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type = (GhostType)gt; Mesh::type_iterator it = mesh.firstType(_all_dimensions, ghost_type, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(_all_dimensions, ghost_type, _ek_not_defined); for (; it != end; ++it) { ElementType type = *it; Array & connectivity_array = mesh.getConnectivity(type, ghost_type); UInt nb_nodes_per_element = connectivity_array.getNbComponent(); Array::vector_iterator conn_it = connectivity_array.begin(nb_nodes_per_element); Array::vector_iterator conn_end = connectivity_array.end(nb_nodes_per_element); for (; conn_it != conn_end; ++conn_it) { for (UInt n = 0; n < nb_nodes_per_element; ++n) { (*conn_it)(n) = new_numbering((*conn_it)(n)); } } } } } this->mesh.sendEvent(remove_nodes); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MeshIgfemSphericalGrowingGel::buildIGFEMMesh() { AKANTU_DEBUG_IN(); NewIGFEMElementsEvent new_elements_event; UInt total_new_elements = 0; Array & new_elements_list = new_elements_event.getList(); Array & old_elements_list = new_elements_event.getOldElementsList(); RemovedElementsEvent removed_elements_event(this->mesh); ChangedElementsEvent changed_elements_event(this->mesh); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; mesh.initElementTypeMapArray(removed_elements_event.getNewNumbering(), 1, dim, ghost_type, false, _ek_not_defined); mesh.initElementTypeMapArray(changed_elements_event.getNewNumbering(), 1, dim, ghost_type, false, _ek_not_defined); /// loop over all types in the mesh for a given ghost type Mesh::type_iterator iit = mesh.firstType(dim, ghost_type, _ek_not_defined); Mesh::type_iterator iend = mesh.lastType(dim, ghost_type, _ek_not_defined); for (; iit != iend; ++iit) { ElementType type = *iit; MeshAbstractIntersector & intersector = *intersectors(type, ghost_type); const Array & new_node_per_elem = intersector.getNewNodePerElem(); UInt n_new_el = 0; Array & connectivity = this->mesh.getConnectivity(type, ghost_type); /// get the connectivities of all types that that may transform Array & connec_igfem_tri_4 = this->mesh.getConnectivity(_igfem_triangle_4, ghost_type); Array & connec_igfem_tri_5 = this->mesh.getConnectivity(_igfem_triangle_5, ghost_type); Array & connec_tri_3 = this->mesh.getConnectivity(_triangle_3, ghost_type); /// create elements to store the newly generated elements Element el_tri_3(_triangle_3, 0, ghost_type, _ek_regular); Element el_igfem_tri_4(_igfem_triangle_4, 0, ghost_type, _ek_igfem); Element el_igfem_tri5(_igfem_triangle_5, 0, ghost_type, _ek_igfem); Array & new_numbering = removed_elements_event.getNewNumbering(type, ghost_type); new_numbering.resize(connectivity.getSize()); /// container for element to be removed Element old_element(type, 0, ghost_type, Mesh::getKind(type)); for (UInt nel = 0; nel != new_node_per_elem.getSize(); ++nel) { /// a former IGFEM triangle is transformed into a regular triangle if ((type != _triangle_3) && (new_node_per_elem(nel, 0) == 0)) { Vector connec_new_elem(3); connec_new_elem(0) = connectivity(nel, 0); connec_new_elem(1) = connectivity(nel, 1); connec_new_elem(2) = connectivity(nel, 2); /// add the new element in the mesh UInt nb_triangles_3 = connec_tri_3.getSize(); connec_tri_3.push_back(connec_new_elem); el_tri_3.element = nb_triangles_3; new_elements_list.push_back(el_tri_3); /// the number of the old element in the mesh old_element.element = nel; old_elements_list.push_back(old_element); new_numbering(nel) = UInt(-1); ++n_new_el; } /// element is enriched else if (new_node_per_elem(nel, 0) != 0) { switch (new_node_per_elem(nel, 0)) { /// new element is of type igfem_triangle_4 case 1: { Vector connec_new_elem(4); switch (new_node_per_elem(nel, 2)) { case 0: connec_new_elem(0) = connectivity(nel, 2); connec_new_elem(1) = connectivity(nel, 0); connec_new_elem(2) = connectivity(nel, 1); break; case 1: connec_new_elem(0) = connectivity(nel, 0); connec_new_elem(1) = connectivity(nel, 1); connec_new_elem(2) = connectivity(nel, 2); break; case 2: connec_new_elem(0) = connectivity(nel, 1); connec_new_elem(1) = connectivity(nel, 2); connec_new_elem(2) = connectivity(nel, 0); break; default: - AKANTU_DEBUG_ERROR("A triangle has only 3 segments not " + AKANTU_ERROR("A triangle has only 3 segments not " << new_node_per_elem(nel, 2)); break; } connec_new_elem(3) = new_node_per_elem(nel, 1); UInt nb_igfem_triangles_4 = connec_igfem_tri_4.getSize(); connec_igfem_tri_4.push_back(connec_new_elem); el_igfem_tri_4.element = nb_igfem_triangles_4; new_elements_list.push_back(el_igfem_tri_4); break; } /// new element is of type igfem_triangle_5 case 2: { Vector connec_new_elem(5); if ((new_node_per_elem(nel, 2) == 0) && (new_node_per_elem(nel, 4) == 1)) { connec_new_elem(0) = connectivity(nel, 1); connec_new_elem(1) = connectivity(nel, 2); connec_new_elem(2) = connectivity(nel, 0); connec_new_elem(3) = new_node_per_elem(nel, 3); connec_new_elem(4) = new_node_per_elem(nel, 1); } else if ((new_node_per_elem(nel, 2) == 0) && (new_node_per_elem(nel, 4) == 2)) { connec_new_elem(0) = connectivity(nel, 0); connec_new_elem(1) = connectivity(nel, 1); connec_new_elem(2) = connectivity(nel, 2); connec_new_elem(3) = new_node_per_elem(nel, 1); connec_new_elem(4) = new_node_per_elem(nel, 3); } else if ((new_node_per_elem(nel, 2) == 1) && (new_node_per_elem(nel, 4) == 2)) { connec_new_elem(0) = connectivity(nel, 2); connec_new_elem(1) = connectivity(nel, 0); connec_new_elem(2) = connectivity(nel, 1); connec_new_elem(3) = new_node_per_elem(nel, 3); connec_new_elem(4) = new_node_per_elem(nel, 1); } else if ((new_node_per_elem(nel, 2) == 1) && (new_node_per_elem(nel, 4) == 0)) { connec_new_elem(0) = connectivity(nel, 1); connec_new_elem(1) = connectivity(nel, 2); connec_new_elem(2) = connectivity(nel, 0); connec_new_elem(3) = new_node_per_elem(nel, 1); connec_new_elem(4) = new_node_per_elem(nel, 3); } else if ((new_node_per_elem(nel, 2) == 2) && (new_node_per_elem(nel, 4) == 0)) { connec_new_elem(0) = connectivity(nel, 0); connec_new_elem(1) = connectivity(nel, 1); connec_new_elem(2) = connectivity(nel, 2); connec_new_elem(3) = new_node_per_elem(nel, 3); connec_new_elem(4) = new_node_per_elem(nel, 1); } else if ((new_node_per_elem(nel, 2) == 2) && (new_node_per_elem(nel, 4) == 1)) { connec_new_elem(0) = connectivity(nel, 2); connec_new_elem(1) = connectivity(nel, 0); connec_new_elem(2) = connectivity(nel, 1); connec_new_elem(3) = new_node_per_elem(nel, 1); connec_new_elem(4) = new_node_per_elem(nel, 3); } else - AKANTU_DEBUG_ERROR("A triangle has only 3 segments (0 to 2) not " + AKANTU_ERROR("A triangle has only 3 segments (0 to 2) not " << new_node_per_elem(nel, 2) << "and" << new_node_per_elem(nel, 4)); UInt nb_igfem_triangles_5 = connec_igfem_tri_5.getSize(); connec_igfem_tri_5.push_back(connec_new_elem); el_igfem_tri5.element = nb_igfem_triangles_5; new_elements_list.push_back(el_igfem_tri5); break; } default: - AKANTU_DEBUG_ERROR("IGFEM cannot add " << new_node_per_elem(nel, 0) + AKANTU_ERROR("IGFEM cannot add " << new_node_per_elem(nel, 0) << " nodes to triangles"); break; } old_element.element = nel; old_elements_list.push_back(old_element); new_numbering(nel) = UInt(-1); ++n_new_el; } } total_new_elements += n_new_el; } } /// update new numbering for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { GhostType ghost_type = *gt; /// loop over all types in the mesh for a given ghost type Mesh::type_iterator iit = mesh.firstType(dim, ghost_type, _ek_not_defined); Mesh::type_iterator iend = mesh.lastType(dim, ghost_type, _ek_not_defined); for (; iit != iend; ++iit) { ElementType type = *iit; Array & new_numbering = removed_elements_event.getNewNumbering(type, ghost_type); UInt el_index = 0; UInt nb_element = this->mesh.getNbElement(type, ghost_type); new_numbering.resize(nb_element); for (UInt e = 0; e < nb_element; ++e) { if (new_numbering(e) != UInt(-1)) { new_numbering(e) = el_index; ++el_index; } } changed_elements_event.getNewNumbering(type, ghost_type) .copy(new_numbering); } } StaticCommunicator::getStaticCommunicator().allReduce(&total_new_elements, 1, _so_sum); if (total_new_elements > 0) { changed_elements_event.getListOld().copy( new_elements_event.getOldElementsList()); changed_elements_event.getListNew().copy(new_elements_event.getList()); this->mesh.sendEvent(changed_elements_event); this->mesh.sendEvent(new_elements_event); Array & removed_list = removed_elements_event.getList(); removed_list.copy(new_elements_event.getOldElementsList()); this->mesh.sendEvent(removed_elements_event); } removeAdditionalNodes(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MeshIgfemSphericalGrowingGel::buildSegmentConnectivityToNodeType() { Mesh mesh_facets(mesh.initMeshFacets()); MeshUtils::buildSegmentToNodeType(mesh, mesh_facets, synchronizer); // only the ghost elements are considered for (UInt g = _not_ghost; g <= _ghost; ++g) { GhostType ghost_type = (GhostType)g; Mesh::type_iterator it = mesh_facets.firstType(1, ghost_type); Mesh::type_iterator end = mesh_facets.lastType(1, ghost_type); for (; it != end; ++it) { ElementType type = *it; const Array & segment_to_nodetype = mesh_facets.getData("segment_to_nodetype", type, ghost_type); const Array & segment_connectivity = mesh_facets.getConnectivity(type, ghost_type); // looping over all the segments Array::const_vector_iterator conn_it = segment_connectivity.begin(segment_connectivity.getNbComponent()); Array::const_vector_iterator conn_end = segment_connectivity.end(segment_connectivity.getNbComponent()); UInt seg_index = 0; UInt n[2]; for (; conn_it != conn_end; ++conn_it, ++seg_index) { Int seg_type = segment_to_nodetype(seg_index); n[0] = (*conn_it)(0); n[1] = (*conn_it)(1); if ((mesh.isMasterNode(n[0]) || mesh.isSlaveNode(n[0])) && (mesh.isMasterNode(n[1]) || mesh.isSlaveNode(n[1]))) { std::sort(n, n + 2); segment_conn_to_node_type[std::make_pair(n[0], n[1])] = seg_type; } } } } } /* -------------------------------------------------------------------------- */ template void MeshIgfemSphericalGrowingGel::updateNodeType( const Array & nodes_list, const Array & new_node_per_elem, ElementType type, GhostType ghost_type) { Array & nodes_type = mesh.getNodesType(); UInt old_nodes = nodes_type.getSize(); UInt new_nodes = nodes_list.getSize(); // exit this function if the simulation in run in serial if (old_nodes == 0 || new_nodes == 0) return; nodes_type.resize(old_nodes + new_nodes); Array checked_node(new_nodes, 1, false); UInt nb_prim_by_el = 0; if ((type == _triangle_3) || (type == _igfem_triangle_4) || (type == _igfem_triangle_5)) { nb_prim_by_el = 3; } else { - AKANTU_DEBUG_ERROR("Not ready for mesh type " << type); + AKANTU_ERROR("Not ready for mesh type " << type); } // determine the node type for the local, master and slave nodes const Array & connectivity = mesh.getConnectivity(type, ghost_type); unordered_map, Int>::type::iterator seg_type_it; unordered_map, Int>::type::iterator seg_type_end = segment_conn_to_node_type.end(); for (UInt el = 0; el < new_node_per_elem.getSize(); ++el) { UInt nb_nodes = new_node_per_elem(el, 0); for (UInt n = 0; n < nb_nodes; ++n) { UInt node_index = new_node_per_elem(el, 2 * n + 1); if (node_index < old_nodes || checked_node(node_index - old_nodes)) continue; // get the elemental connectivity of the segment associated to the node UInt segment_index = new_node_per_elem(el, 2 * n + 2); UInt extreme_nodes[2]; extreme_nodes[0] = segment_index; extreme_nodes[1] = segment_index + 1; if (extreme_nodes[1] == nb_prim_by_el) extreme_nodes[1] = 0; // get the connectivity of the segment extreme_nodes[0] = connectivity(el, extreme_nodes[0]); extreme_nodes[1] = connectivity(el, extreme_nodes[1]); // if one extreme nodes is pure ghost, then also the new node is pure // ghost if (mesh.isPureGhostNode(extreme_nodes[0]) || mesh.isPureGhostNode(extreme_nodes[1])) nodes_type(node_index) = -3; // if on of the two extreme nodes is local, then also the new node is // local else if (mesh.isLocalNode(extreme_nodes[0]) || mesh.isLocalNode(extreme_nodes[1])) nodes_type(node_index) = -1; // otherwise use the value stored in the map else { std::sort(extreme_nodes, extreme_nodes + 2); seg_type_it = segment_conn_to_node_type.find( std::make_pair(extreme_nodes[0], extreme_nodes[1])); AKANTU_DEBUG_ASSERT(seg_type_it != seg_type_end, "Segment not found"); nodes_type(node_index) = seg_type_it->second; } checked_node(node_index - old_nodes) = true; } } AKANTU_DEBUG_ASSERT(std::accumulate(checked_node.begin(), checked_node.end(), 0) == Int(checked_node.getSize()), "Not all new nodes were assigned a node type"); } __END_AKANTU__ diff --git a/extra_packages/igfem/test/test_fe_engine/igfem_mesh_generation.cc b/extra_packages/igfem/test/test_fe_engine/igfem_mesh_generation.cc index d5e1b3a6c..32770de4a 100644 --- a/extra_packages/igfem/test/test_fe_engine/igfem_mesh_generation.cc +++ b/extra_packages/igfem/test/test_fe_engine/igfem_mesh_generation.cc @@ -1,115 +1,115 @@ /** * @file igfem_mesh_generation.cc * @author Aurelia Isabel Cuba Ramos * @date Fri Oct 30 14:55:51 2015 * * @brief function to generate a IGFEM mesh for fe_engine tests * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "mesh.hh" #include /* -------------------------------------------------------------------------- */ using namespace akantu; void generateIGFEMMesh(const ElementType type, Mesh & mesh, const std::string & filename) { std::ifstream infile; infile.open(filename.c_str()); if (!infile.good()) { - AKANTU_DEBUG_ERROR("Cannot open file " << filename); + AKANTU_ERROR("Cannot open file " << filename); } UInt current_line = 0; std::string line; UInt first_node_number = std::numeric_limits::max(); UInt last_node_number = 0; while (infile.good()) { std::getline(infile, line); current_line++; /// read all nodes if (line == "$Nodes" || line == "$NOD") { UInt nb_nodes; std::getline(infile, line); std::stringstream sstr(line); sstr >> nb_nodes; current_line++; Array & nodes = const_cast &>(mesh.getNodes()); nodes.resize(nb_nodes); UInt index; Real coord[3]; UInt spatial_dimension = nodes.getNbComponent(); /// for each node, read the coordinates for (UInt i = 0; i < nb_nodes; ++i) { UInt offset = i * spatial_dimension; std::getline(infile, line); std::stringstream sstr_node(line); sstr_node >> index >> coord[0] >> coord[1] >> coord[2]; current_line++; first_node_number = std::min(first_node_number, index); last_node_number = std::max(last_node_number, index); /// read the coordinates for (UInt j = 0; j < spatial_dimension; ++j) nodes.storage()[offset + j] = coord[j]; } std::getline(infile, line); /// the end of block line } /// read all elements if (line == "$Elements" || line == "$ELM") { mesh.addConnectivityType(type); Array & connectivity = const_cast &>(mesh.getConnectivity(type)); UInt node_per_element = connectivity.getNbComponent(); UInt nb_elements = 0; std::getline(infile, line); std::stringstream sstr(line); sstr >> nb_elements; current_line++; for (UInt i = 0; i < nb_elements; ++i) { std::getline(infile, line); std::stringstream sstr_elem(line); current_line++; Vector local_connect(node_per_element); for (UInt j = 0; j < node_per_element; ++j) { UInt node_index; sstr_elem >> node_index; AKANTU_DEBUG_ASSERT(node_index <= last_node_number, "Node number not in range : line " << current_line); node_index -= first_node_number; local_connect(j) = node_index; } connectivity.push_back(local_connect); } std::getline(infile, line); /// the end of block line } } } diff --git a/extra_packages/traction-at-split-node-contact/src/boundary_conditions/force_based_dirichlet.hh b/extra_packages/traction-at-split-node-contact/src/boundary_conditions/force_based_dirichlet.hh index e57d1a007..014e6959f 100644 --- a/extra_packages/traction-at-split-node-contact/src/boundary_conditions/force_based_dirichlet.hh +++ b/extra_packages/traction-at-split-node-contact/src/boundary_conditions/force_based_dirichlet.hh @@ -1,132 +1,128 @@ /** * @file force_based_dirichlet.hh * * @author Dana Christen * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief dirichlet boundary condition that tries * to keep the force at a given value * * @section LICENSE * * Copyright (©) 2015-2018 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 __AST_FORCE_BASED_DIRICHLET_HH__ #define __AST_FORCE_BASED_DIRICHLET_HH__ // akantu #include "aka_common.hh" namespace akantu { /* -------------------------------------------------------------------------- */ class ForceBasedDirichlet : public BC::Dirichlet::IncrementValue { protected: typedef const Array * RealArrayPtr; typedef const Array * IntArrayPtr; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: ForceBasedDirichlet(SolidMechanicsModel & model, BC::Axis ax, Real target_f, Real mass = 0.) : IncrementValue(0., ax), model(model), mass(mass), velocity(0.), target_force(target_f), total_residual(0.) {} virtual ~ForceBasedDirichlet() {} /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void updateTotalResidual() { - SubBoundarySet::iterator it = this->subboundaries.begin(); - SubBoundarySet::iterator end = this->subboundaries.end(); this->total_residual = 0.; - for (; it != end; ++it) { - this->total_residual += integrateResidual(*it, this->model, this->axis); + for (auto && subboundary : this->subboundaries) { + this->total_residual += integrateResidual(subboundary, this->model, this->axis); } } virtual Real update() { AKANTU_DEBUG_IN(); this->updateTotalResidual(); Real total_force = this->target_force + this->total_residual; Real a = total_force / this->mass; Real dt = model.getTimeStep(); this->velocity += 0.5 * dt * a; this->value = this->velocity * dt + 0.5 * dt * dt * a; // increment position dx this->velocity += 0.5 * dt * a; AKANTU_DEBUG_OUT(); return this->total_residual; } Real applyYourself() { AKANTU_DEBUG_IN(); Real reaction = this->update(); - SubBoundarySet::iterator it = this->subboundaries.begin(); - SubBoundarySet::iterator end = this->subboundaries.end(); - for (; it != end; ++it) { - this->model.applyBC(*this, *it); + for (auto && subboundary : this->subboundaries) { + this->model.applyBC(*this, subboundary); } AKANTU_DEBUG_OUT(); return reaction; } /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_SET_MACRO(Mass, mass, Real); AKANTU_SET_MACRO(TargetForce, target_force, Real); void insertSubBoundary(const std::string & sb_name) { this->subboundaries.insert(sb_name); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ typedef std::set SubBoundarySet; protected: SolidMechanicsModel & model; SubBoundarySet subboundaries; Real mass; Real velocity; Real target_force; Real total_residual; }; } // namespace akantu #endif /* __AST_FORCE_BASED_DIRICHLET_HH__ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_coulomb.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_coulomb.hh index 0cc6df6dd..634efba72 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_coulomb.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_coulomb.hh @@ -1,108 +1,108 @@ /** * @file ntn_friclaw_coulomb.hh * * @author David Simon Kammer * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Feb 23 2018 * * @brief coulomb friction with \mu_s = \mu_k (constant) * * @section LICENSE * * Copyright (©) 2015-2018 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 __AST_NTN_FRICLAW_COULOMB_HH__ #define __AST_NTN_FRICLAW_COULOMB_HH__ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_fricreg_no_regularisation.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template class NTNFricLawCoulomb : public Regularisation { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - NTNFricLawCoulomb(NTNBaseContact * contact, const FrictionID & id = "coulomb", + NTNFricLawCoulomb(NTNBaseContact & contact, const ID & id = "coulomb", const MemoryID & memory_id = 0); virtual ~NTNFricLawCoulomb(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// register synchronizedarrays for sync virtual void registerSynchronizedArray(SynchronizedArrayBase & array); /// dump restart file virtual void dumpRestart(const std::string & file_name) const; /// read restart file virtual void readRestart(const std::string & file_name); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; protected: /// compute frictional strength according to friction law virtual void computeFrictionalStrength(); /* ------------------------------------------------------------------------ */ /* Dumpable */ /* ------------------------------------------------------------------------ */ public: virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: // friction coefficient SynchronizedArray mu; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ /// standard output stream operator template inline std::ostream & operator<<(std::ostream & stream, const NTNFricLawCoulomb & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "ntn_friclaw_coulomb_tmpl.hh" #endif /* __AST_NTN_FRICLAW_COULOMB_HH__ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_coulomb_tmpl.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_coulomb_tmpl.hh index d012fb969..112f9fe4c 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_coulomb_tmpl.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_coulomb_tmpl.hh @@ -1,168 +1,168 @@ /** * @file ntn_friclaw_coulomb_tmpl.hh * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief implementation of coulomb friction * * @section LICENSE * * Copyright (©) 2015-2018 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 "dumper_nodal_field.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template -NTNFricLawCoulomb::NTNFricLawCoulomb(NTNBaseContact * contact, - const FrictionID & id, +NTNFricLawCoulomb::NTNFricLawCoulomb(NTNBaseContact & contact, + const ID & id, const MemoryID & memory_id) : Regularisation(contact, id, memory_id), mu(0, 1, 0., id + ":mu", 0., "mu") { AKANTU_DEBUG_IN(); Regularisation::registerSynchronizedArray(this->mu); this->registerParam("mu", this->mu, _pat_parsmod, "friction coefficient"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawCoulomb::computeFrictionalStrength() { AKANTU_DEBUG_IN(); // get contact arrays const SynchronizedArray & is_in_contact = this->internalGetIsInContact(); const SynchronizedArray & pressure = this->internalGetContactPressure(); // array to fill SynchronizedArray & strength = this->internalGetFrictionalStrength(); - UInt nb_contact_nodes = this->contact->getNbContactNodes(); + UInt nb_contact_nodes = this->contact.getNbContactNodes(); for (UInt n = 0; n < nb_contact_nodes; ++n) { // node pair is NOT in contact if (!is_in_contact(n)) strength(n) = 0.; // node pair is in contact else { // compute frictional strength strength(n) = this->mu(n) * pressure(n); } } Regularisation::computeFrictionalStrength(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawCoulomb::registerSynchronizedArray( SynchronizedArrayBase & array) { AKANTU_DEBUG_IN(); this->mu.registerDependingArray(array); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawCoulomb::dumpRestart( const std::string & file_name) const { AKANTU_DEBUG_IN(); this->mu.dumpRestartFile(file_name); Regularisation::dumpRestart(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawCoulomb::readRestart( const std::string & file_name) { AKANTU_DEBUG_IN(); this->mu.readRestartFile(file_name); Regularisation::readRestart(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawCoulomb::printself(std::ostream & stream, int indent) const { AKANTU_DEBUG_IN(); std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "NTNFricLawCoulomb [" << std::endl; Regularisation::printself(stream, ++indent); stream << space << "]" << std::endl; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawCoulomb::addDumpFieldToDumper( const std::string & dumper_name, const std::string & field_id) { AKANTU_DEBUG_IN(); #ifdef AKANTU_USE_IOHELPER // const SynchronizedArray * nodal_filter = - // &(this->contact->getSlaves()); + // &(this->contact.getSlaves()); if (field_id == "mu") { this->internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->mu.getArray())); } /* else if (field_id == "frictional_contact_pressure") { this->internalAddDumpFieldToDumper(dumper_name, field_id, new DumperIOHelper::NodalField(this->frictional_contact_pressure.getArray())); } */ else { Regularisation::addDumpFieldToDumper(dumper_name, field_id); } #endif AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_cohesive.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_cohesive.hh index a6bb69a06..5a29f3e87 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_cohesive.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_cohesive.hh @@ -1,115 +1,115 @@ /** * @file ntn_friclaw_linear_cohesive.hh * * @author David Simon Kammer * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Feb 23 2018 * * @brief linear cohesive law * * @section LICENSE * * Copyright (©) 2015-2018 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 __AST_NTN_FRICLAW_LINEAR_COHESIVE_HH__ #define __AST_NTN_FRICLAW_LINEAR_COHESIVE_HH__ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_fricreg_no_regularisation.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template class NTNFricLawLinearCohesive : public Regularisation { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - NTNFricLawLinearCohesive(NTNBaseContact * contact, - const FrictionID & id = "linear_cohesive", + NTNFricLawLinearCohesive(NTNBaseContact & contact, + const ID & id = "linear_cohesive", const MemoryID & memory_id = 0); virtual ~NTNFricLawLinearCohesive(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// register synchronizedarrays for sync virtual void registerSynchronizedArray(SynchronizedArrayBase & array); /// dump restart file virtual void dumpRestart(const std::string & file_name) const; /// read restart file virtual void readRestart(const std::string & file_name); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; protected: /// compute frictional strength according to friction law virtual void computeFrictionalStrength(); /* ------------------------------------------------------------------------ */ /* Dumpable */ /* ------------------------------------------------------------------------ */ public: virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: // fracture energy SynchronizedArray G_c; // peak value of cohesive law SynchronizedArray tau_c; // residual value of cohesive law (for slip > d_c) SynchronizedArray tau_r; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ /// standard output stream operator template inline std::ostream & operator<<(std::ostream & stream, const NTNFricLawLinearCohesive & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "ntn_friclaw_linear_cohesive_tmpl.hh" #endif /* __AST_NTN_FRICLAW_LINEAR_COHESIVE_HH__ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_cohesive_tmpl.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_cohesive_tmpl.hh index bd9fe4f2d..724a9cd79 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_cohesive_tmpl.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_cohesive_tmpl.hh @@ -1,189 +1,189 @@ /** * @file ntn_friclaw_linear_cohesive_tmpl.hh * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief implementation of linear cohesive law * * @section LICENSE * * Copyright (©) 2015-2018 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 "dumper_text.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template NTNFricLawLinearCohesive::NTNFricLawLinearCohesive( - NTNBaseContact * contact, const FrictionID & id, const MemoryID & memory_id) + NTNBaseContact & contact, const ID & id, const MemoryID & memory_id) : Regularisation(contact, id, memory_id), G_c(0, 1, 0., id + ":G_c", 0., "G_c"), tau_c(0, 1, 0., id + ":tau_c", 0., "tau_c"), tau_r(0, 1, 0., id + ":tau_r", 0., "tau_r") { AKANTU_DEBUG_IN(); Regularisation::registerSynchronizedArray(this->G_c); Regularisation::registerSynchronizedArray(this->tau_c); Regularisation::registerSynchronizedArray(this->tau_r); this->registerParam("G_c", this->G_c, _pat_parsmod, "fracture energy"); this->registerParam("tau_c", this->tau_c, _pat_parsmod, "peak shear strength"); this->registerParam("tau_r", this->tau_r, _pat_parsmod, "residual shear strength"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearCohesive::computeFrictionalStrength() { AKANTU_DEBUG_IN(); // get arrays const SynchronizedArray & is_in_contact = this->internalGetIsInContact(); // const SynchronizedArray & slip = this->internalGetSlip(); const SynchronizedArray & slip = this->internalGetCumulativeSlip(); // array to fill SynchronizedArray & strength = this->internalGetFrictionalStrength(); - UInt nb_contact_nodes = this->contact->getNbContactNodes(); + UInt nb_contact_nodes = this->contact.getNbContactNodes(); for (UInt n = 0; n < nb_contact_nodes; ++n) { // node pair is NOT in contact if (!is_in_contact(n)) strength(n) = 0.; // node pair is in contact else { if (this->G_c(n) == 0.) { // strength(n) = 0.; strength(n) = this->tau_r(n); } else { Real slope = (this->tau_c(n) - this->tau_r(n)) * (this->tau_c(n) - this->tau_r(n)) / (2 * this->G_c(n)); // no strength < tau_r strength(n) = std::max(this->tau_c(n) - slope * slip(n), this->tau_r(n)); // strength(n) = std::max(this->tau_c(n) - slope * slip(n), 0.); // no // negative strength } } } Regularisation::computeFrictionalStrength(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearCohesive::registerSynchronizedArray( SynchronizedArrayBase & array) { AKANTU_DEBUG_IN(); this->G_c.registerDependingArray(array); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearCohesive::dumpRestart( const std::string & file_name) const { AKANTU_DEBUG_IN(); this->G_c.dumpRestartFile(file_name); this->tau_c.dumpRestartFile(file_name); this->tau_r.dumpRestartFile(file_name); Regularisation::dumpRestart(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearCohesive::readRestart( const std::string & file_name) { AKANTU_DEBUG_IN(); this->G_c.readRestartFile(file_name); this->tau_c.readRestartFile(file_name); this->tau_r.readRestartFile(file_name); Regularisation::readRestart(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearCohesive::printself(std::ostream & stream, int indent) const { AKANTU_DEBUG_IN(); std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "NTNFricLawLinearCohesive [" << std::endl; Regularisation::printself(stream, ++indent); stream << space << "]" << std::endl; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearCohesive::addDumpFieldToDumper( const std::string & dumper_name, const std::string & field_id) { AKANTU_DEBUG_IN(); #ifdef AKANTU_USE_IOHELPER // const SynchronizedArray * nodal_filter = - // &(this->contact->getSlaves()); + // &(this->contact.getSlaves()); if (field_id == "G_c") { this->internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->G_c.getArray())); } else if (field_id == "tau_c") { this->internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->tau_c.getArray())); } else if (field_id == "tau_r") { this->internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->tau_r.getArray())); } else { Regularisation::addDumpFieldToDumper(dumper_name, field_id); } #endif AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening.hh index b95d989c0..347a9dbdd 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening.hh @@ -1,117 +1,117 @@ /** * @file ntn_friclaw_linear_slip_weakening.hh * * @author David Simon Kammer * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Feb 23 2018 * * @brief linear slip weakening * * @section LICENSE * * Copyright (©) 2015-2018 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 __AST_NTN_FRICLAW_LINEAR_SLIP_WEAKENING_HH__ #define __AST_NTN_FRICLAW_LINEAR_SLIP_WEAKENING_HH__ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_friclaw_coulomb.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template class NTNFricLawLinearSlipWeakening : public NTNFricLawCoulomb { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - NTNFricLawLinearSlipWeakening(NTNBaseContact * contact, - const FrictionID & id = "linear_slip_weakening", + NTNFricLawLinearSlipWeakening(NTNBaseContact & contact, + const ID & id = "linear_slip_weakening", const MemoryID & memory_id = 0); virtual ~NTNFricLawLinearSlipWeakening(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// register synchronizedarrays for sync virtual void registerSynchronizedArray(SynchronizedArrayBase & array); /// dump restart file virtual void dumpRestart(const std::string & file_name) const; /// read restart file virtual void readRestart(const std::string & file_name); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; protected: /// compute frictional strength according to friction law virtual void computeFrictionalStrength(); /// computes the friction coefficient as a function of slip virtual void computeFrictionCoefficient(); /* ------------------------------------------------------------------------ */ /* Dumpable */ /* ------------------------------------------------------------------------ */ public: virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: // static coefficient of friction SynchronizedArray mu_s; // kinetic coefficient of friction SynchronizedArray mu_k; // Dc the length over which slip weakening happens SynchronizedArray d_c; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ /// standard output stream operator template inline std::ostream & operator<<(std::ostream & stream, const NTNFricLawLinearSlipWeakening & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "ntn_friclaw_linear_slip_weakening_tmpl.hh" #endif /* __AST_NTN_FRICLAW_LINEAR_SLIP_WEAKENING_HH__ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening_no_healing.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening_no_healing.hh index 9711597b5..14f8d1de4 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening_no_healing.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening_no_healing.hh @@ -1,96 +1,96 @@ /** * @file ntn_friclaw_linear_slip_weakening_no_healing.hh * * @author David Simon Kammer * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Feb 23 2018 * * @brief linear slip weakening * * @section LICENSE * * Copyright (©) 2015-2018 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 __AST_NTN_FRICLAW_LINEAR_SLIP_WEAKENING_NO_HEALING_HH__ #define __AST_NTN_FRICLAW_LINEAR_SLIP_WEAKENING_NO_HEALING_HH__ /* -------------------------------------------------------------------------- */ #include "ntn_friclaw_linear_slip_weakening.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template class NTNFricLawLinearSlipWeakeningNoHealing : public NTNFricLawLinearSlipWeakening { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: NTNFricLawLinearSlipWeakeningNoHealing( - NTNBaseContact * contact, - const FrictionID & id = "linear_slip_weakening_no_healing", + NTNBaseContact & contact, + const ID & id = "linear_slip_weakening_no_healing", const MemoryID & memory_id = 0); virtual ~NTNFricLawLinearSlipWeakeningNoHealing(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; protected: /// computes the friction coefficient as a function of slip virtual void computeFrictionCoefficient(); /* ------------------------------------------------------------------------ */ /* Dumpable */ /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ /// standard output stream operator template inline std::ostream & operator<<( std::ostream & stream, const NTNFricLawLinearSlipWeakeningNoHealing & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "ntn_friclaw_linear_slip_weakening_no_healing_tmpl.hh" #endif /* __AST_NTN_FRICLAW_LINEAR_SLIP_WEAKENING_NO_HEALING_HH__ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening_no_healing_tmpl.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening_no_healing_tmpl.hh index a58a7d954..2341d348b 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening_no_healing_tmpl.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening_no_healing_tmpl.hh @@ -1,87 +1,87 @@ /** * @file ntn_friclaw_linear_slip_weakening_no_healing_tmpl.hh * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief implementation of linear slip weakening * * @section LICENSE * * Copyright (©) 2015-2018 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 . * */ /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template NTNFricLawLinearSlipWeakeningNoHealing:: - NTNFricLawLinearSlipWeakeningNoHealing(NTNBaseContact * contact, - const FrictionID & id, + NTNFricLawLinearSlipWeakeningNoHealing(NTNBaseContact & contact, + const ID & id, const MemoryID & memory_id) : NTNFricLawLinearSlipWeakening(contact, id, memory_id) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearSlipWeakeningNoHealing< Regularisation>::computeFrictionCoefficient() { AKANTU_DEBUG_IN(); // get arrays const SynchronizedArray & slip = this->internalGetCumulativeSlip(); - UInt nb_contact_nodes = this->contact->getNbContactNodes(); + UInt nb_contact_nodes = this->contact.getNbContactNodes(); for (UInt n = 0; n < nb_contact_nodes; ++n) { if (slip(n) >= this->d_c(n)) { this->mu(n) = this->mu_k(n); } else { // mu = mu_k + (1 - slip / Dc) * (mu_s - mu_k) this->mu(n) = this->mu_k(n) + (1 - (slip(n) / this->d_c(n))) * (this->mu_s(n) - this->mu_k(n)); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearSlipWeakeningNoHealing::printself( std::ostream & stream, int indent) const { AKANTU_DEBUG_IN(); std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "NTNFricLawLinearSlipWeakeningNoHealing [" << std::endl; NTNFricLawLinearSlipWeakening::printself(stream, ++indent); stream << space << "]" << std::endl; AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening_tmpl.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening_tmpl.hh index f54a6c269..c92037880 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening_tmpl.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_laws/ntn_friclaw_linear_slip_weakening_tmpl.hh @@ -1,191 +1,191 @@ /** * @file ntn_friclaw_linear_slip_weakening_tmpl.hh * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief implementation of linear slip weakening * * @section LICENSE * * Copyright (©) 2015-2018 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 "dumper_text.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template NTNFricLawLinearSlipWeakening::NTNFricLawLinearSlipWeakening( - NTNBaseContact * contact, const FrictionID & id, const MemoryID & memory_id) + NTNBaseContact & contact, const ID & id, const MemoryID & memory_id) : NTNFricLawCoulomb(contact, id, memory_id), mu_s(0, 1, 0., id + ":mu_s", 0., "mu_s"), mu_k(0, 1, 0., id + ":mu_k", 0., "mu_k"), d_c(0, 1, 0., id + ":d_c", 0., "d_c") { AKANTU_DEBUG_IN(); NTNFricLawCoulomb::registerSynchronizedArray(this->mu_s); NTNFricLawCoulomb::registerSynchronizedArray(this->mu_k); NTNFricLawCoulomb::registerSynchronizedArray(this->d_c); this->registerParam("mu_s", this->mu_s, _pat_parsmod, "static friction coefficient"); this->registerParam("mu_k", this->mu_k, _pat_parsmod, "kinetic friction coefficient"); this->registerParam("d_c", this->d_c, _pat_parsmod, "slip weakening length"); - this->setParamAccessType("mu", _pat_readable); + this->setParameterAccessType("mu", _pat_readable); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearSlipWeakening< Regularisation>::computeFrictionalStrength() { AKANTU_DEBUG_IN(); computeFrictionCoefficient(); NTNFricLawCoulomb::computeFrictionalStrength(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearSlipWeakening< Regularisation>::computeFrictionCoefficient() { AKANTU_DEBUG_IN(); // get arrays const SynchronizedArray & stick = this->internalGetIsSticking(); const SynchronizedArray & slip = this->internalGetSlip(); - UInt nb_contact_nodes = this->contact->getNbContactNodes(); + UInt nb_contact_nodes = this->contact.getNbContactNodes(); for (UInt n = 0; n < nb_contact_nodes; ++n) { if (stick(n)) { this->mu(n) = this->mu_s(n); } else { if (slip(n) >= this->d_c(n)) { this->mu(n) = this->mu_k(n); } else { // mu = mu_k + (1 - slip / Dc) * (mu_s - mu_k) this->mu(n) = this->mu_k(n) + (1 - (slip(n) / this->d_c(n))) * (this->mu_s(n) - this->mu_k(n)); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearSlipWeakening::registerSynchronizedArray( SynchronizedArrayBase & array) { AKANTU_DEBUG_IN(); this->mu_s.registerDependingArray(array); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearSlipWeakening::dumpRestart( const std::string & file_name) const { AKANTU_DEBUG_IN(); this->mu_s.dumpRestartFile(file_name); this->mu_k.dumpRestartFile(file_name); this->d_c.dumpRestartFile(file_name); NTNFricLawCoulomb::dumpRestart(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearSlipWeakening::readRestart( const std::string & file_name) { AKANTU_DEBUG_IN(); this->mu_s.readRestartFile(file_name); this->mu_k.readRestartFile(file_name); this->d_c.readRestartFile(file_name); NTNFricLawCoulomb::readRestart(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearSlipWeakening::printself( std::ostream & stream, int indent) const { AKANTU_DEBUG_IN(); std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "NTNFricLawLinearSlipWeakening [" << std::endl; NTNFricLawCoulomb::printself(stream, ++indent); stream << space << "]" << std::endl; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void NTNFricLawLinearSlipWeakening::addDumpFieldToDumper( const std::string & dumper_name, const std::string & field_id) { AKANTU_DEBUG_IN(); #ifdef AKANTU_USE_IOHELPER // const SynchronizedArray * nodal_filter = - // &(this->contact->getSlaves()); + // &(this->contact.getSlaves()); if (field_id == "mu_s") { this->internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->mu_s.getArray())); } else if (field_id == "mu_k") { this->internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->mu_k.getArray())); } else if (field_id == "d_c") { this->internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->d_c.getArray())); } else { NTNFricLawCoulomb::addDumpFieldToDumper(dumper_name, field_id); } #endif AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_no_regularisation.cc b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_no_regularisation.cc index 60a69dd3b..323e4cae3 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_no_regularisation.cc +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_no_regularisation.cc @@ -1,168 +1,168 @@ /** * @file ntn_fricreg_no_regularisation.cc * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief implementation of no regularisation * * @section LICENSE * * Copyright (©) 2015-2018 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 . * */ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_fricreg_no_regularisation.hh" #include "dumper_nodal_field.hh" #include "dumper_text.hh" namespace akantu { /* -------------------------------------------------------------------------- */ NTNFricRegNoRegularisation::NTNFricRegNoRegularisation( - NTNBaseContact * contact, const FrictionID & id, const MemoryID & memory_id) + NTNBaseContact & contact, const ID & id, const MemoryID & memory_id) : NTNBaseFriction(contact, id, memory_id), frictional_contact_pressure(0, 1, 0., id + ":frictional_contact_pressure", 0., "frictional_contact_pressure") { AKANTU_DEBUG_IN(); NTNBaseFriction::registerSynchronizedArray(this->frictional_contact_pressure); this->registerParam("frictional_contact_pressure", this->frictional_contact_pressure, _pat_internal, "contact pressure used for friction law"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ const SynchronizedArray & NTNFricRegNoRegularisation::internalGetContactPressure() { AKANTU_DEBUG_IN(); this->computeFrictionalContactPressure(); AKANTU_DEBUG_OUT(); return this->frictional_contact_pressure; } /* -------------------------------------------------------------------------- */ void NTNFricRegNoRegularisation::computeFrictionalContactPressure() { AKANTU_DEBUG_IN(); - SolidMechanicsModel & model = this->contact->getModel(); + SolidMechanicsModel & model = this->contact.getModel(); UInt dim = model.getSpatialDimension(); // get contact arrays const SynchronizedArray & is_in_contact = this->internalGetIsInContact(); - const Array & pressure = this->contact->getContactPressure().getArray(); + const Array & pressure = this->contact.getContactPressure().getArray(); Array::const_iterator> it = pressure.begin(dim); - UInt nb_contact_nodes = this->contact->getNbContactNodes(); + UInt nb_contact_nodes = this->contact.getNbContactNodes(); for (UInt n = 0; n < nb_contact_nodes; ++n) { // node pair is NOT in contact if (!is_in_contact(n)) this->frictional_contact_pressure(n) = 0.; // node pair is in contact else { // compute frictional contact pressure const Vector & pres = it[n]; this->frictional_contact_pressure(n) = pres.norm(); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegNoRegularisation::registerSynchronizedArray( SynchronizedArrayBase & array) { AKANTU_DEBUG_IN(); this->frictional_contact_pressure.registerDependingArray(array); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegNoRegularisation::dumpRestart( const std::string & file_name) const { AKANTU_DEBUG_IN(); this->frictional_contact_pressure.dumpRestartFile(file_name); NTNBaseFriction::dumpRestart(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegNoRegularisation::readRestart(const std::string & file_name) { AKANTU_DEBUG_IN(); this->frictional_contact_pressure.readRestartFile(file_name); NTNBaseFriction::readRestart(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegNoRegularisation::printself(std::ostream & stream, int indent) const { AKANTU_DEBUG_IN(); std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "NTNFricRegNoRegularisation [" << std::endl; NTNBaseFriction::printself(stream, ++indent); stream << space << "]" << std::endl; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegNoRegularisation::addDumpFieldToDumper( const std::string & dumper_name, const std::string & field_id) { AKANTU_DEBUG_IN(); #ifdef AKANTU_USE_IOHELPER // const SynchronizedArray * nodal_filter = - // &(this->contact->getSlaves()); + // &(this->contact.getSlaves()); if (field_id == "frictional_contact_pressure") { this->internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField( this->frictional_contact_pressure.getArray())); } else { NTNBaseFriction::addDumpFieldToDumper(dumper_name, field_id); } #endif AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_no_regularisation.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_no_regularisation.hh index 8bf0f4173..9d50e51c5 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_no_regularisation.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_no_regularisation.hh @@ -1,134 +1,134 @@ /** * @file ntn_fricreg_no_regularisation.hh * * @author David Simon Kammer * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Feb 23 2018 * * @brief regularisation that does nothing * * @section LICENSE * * Copyright (©) 2015-2018 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 __AST_NTN_FRICREG_NO_REGULARISATION_HH__ #define __AST_NTN_FRICREG_NO_REGULARISATION_HH__ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_base_friction.hh" namespace akantu { /* -------------------------------------------------------------------------- */ class NTNFricRegNoRegularisation : public NTNBaseFriction { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - NTNFricRegNoRegularisation(NTNBaseContact * contact, - const FrictionID & id = "no_regularisation", + NTNFricRegNoRegularisation(NTNBaseContact & contact, + const ID & id = "no_regularisation", const MemoryID & memory_id = 0); virtual ~NTNFricRegNoRegularisation(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// set to steady state for no regularisation -> do nothing virtual void setToSteadyState(){}; virtual void registerSynchronizedArray(SynchronizedArrayBase & array); virtual void dumpRestart(const std::string & file_name) const; virtual void readRestart(const std::string & file_name); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; protected: virtual void computeFrictionalContactPressure(); /// compute frictional strength according to friction law virtual void computeFrictionalStrength(){}; /* ------------------------------------------------------------------------ */ /* Dumpable */ /* ------------------------------------------------------------------------ */ public: virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: protected: /// get the is_in_contact array virtual const SynchronizedArray & internalGetIsInContact() { - return this->contact->getIsInContact(); + return this->contact.getIsInContact(); }; /// get the contact pressure (the norm: scalar value) virtual const SynchronizedArray & internalGetContactPressure(); /// get the frictional strength array virtual SynchronizedArray & internalGetFrictionalStrength() { return this->frictional_strength; }; /// get the is_sticking array virtual SynchronizedArray & internalGetIsSticking() { return this->is_sticking; } /// get the slip array virtual SynchronizedArray & internalGetSlip() { return this->slip; } /// get the slip array virtual SynchronizedArray & internalGetCumulativeSlip() { return this->cumulative_slip; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: // contact pressure (absolut value) for computation of friction SynchronizedArray frictional_contact_pressure; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "ntn_fricreg_no_regularisation_inline_impl.cc" /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const NTNFricRegNoRegularisation & _this) { _this.printself(stream); return stream; } } // namespace akantu #endif /* __AST_NTN_FRICREG_NO_REGULARISATION_HH__ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_rubin_ampuero.cc b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_rubin_ampuero.cc index 20a82af49..57c2d4ade 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_rubin_ampuero.cc +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_rubin_ampuero.cc @@ -1,176 +1,176 @@ /** * @file ntn_fricreg_rubin_ampuero.cc * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief implementation of no regularisation * * @section LICENSE * * Copyright (©) 2015-2018 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 . * */ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_fricreg_rubin_ampuero.hh" #include "dumper_nodal_field.hh" #include "dumper_text.hh" namespace akantu { /* -------------------------------------------------------------------------- */ -NTNFricRegRubinAmpuero::NTNFricRegRubinAmpuero(NTNBaseContact * contact, - const FrictionID & id, +NTNFricRegRubinAmpuero::NTNFricRegRubinAmpuero(NTNBaseContact & contact, + const ID & id, const MemoryID & memory_id) : NTNFricRegNoRegularisation(contact, id, memory_id), t_star(0, 1, 0., id + ":t_star", 0., "t_star") { AKANTU_DEBUG_IN(); NTNFricRegNoRegularisation::registerSynchronizedArray(this->t_star); this->registerParam("t_star", this->t_star, _pat_parsmod, "time scale of regularization"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ const SynchronizedArray & NTNFricRegRubinAmpuero::internalGetContactPressure() { AKANTU_DEBUG_IN(); - SolidMechanicsModel & model = this->contact->getModel(); + SolidMechanicsModel & model = this->contact.getModel(); UInt dim = model.getSpatialDimension(); Real delta_t = model.getTimeStep(); // get contact arrays const SynchronizedArray & is_in_contact = this->internalGetIsInContact(); - const Array & pressure = this->contact->getContactPressure().getArray(); + const Array & pressure = this->contact.getContactPressure().getArray(); Array::const_iterator> it = pressure.begin(dim); - UInt nb_contact_nodes = this->contact->getNbContactNodes(); + UInt nb_contact_nodes = this->contact.getNbContactNodes(); for (UInt n = 0; n < nb_contact_nodes; ++n) { // node pair is NOT in contact if (!is_in_contact(n)) this->frictional_contact_pressure(n) = 0.; // if t_star is too small compute like Coulomb friction (without // regularization) else if (Math::are_float_equal(this->t_star(n), 0.)) { const Vector & pres = it[n]; this->frictional_contact_pressure(n) = pres.norm(); } else { // compute frictional contact pressure // backward euler method: first order implicit numerical integration // method // \reg_pres_n+1 = (\reg_pres_n + \delta_t / \t_star * \cur_pres) // / (1 + \delta_t / \t_star) Real alpha = delta_t / this->t_star(n); const Vector & pres = it[n]; this->frictional_contact_pressure(n) += alpha * pres.norm(); this->frictional_contact_pressure(n) /= 1 + alpha; } } AKANTU_DEBUG_OUT(); return this->frictional_contact_pressure; } /* -------------------------------------------------------------------------- */ void NTNFricRegRubinAmpuero::setToSteadyState() { AKANTU_DEBUG_IN(); NTNFricRegNoRegularisation::computeFrictionalContactPressure(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegRubinAmpuero::registerSynchronizedArray( SynchronizedArrayBase & array) { AKANTU_DEBUG_IN(); this->t_star.registerDependingArray(array); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegRubinAmpuero::dumpRestart(const std::string & file_name) const { AKANTU_DEBUG_IN(); this->t_star.dumpRestartFile(file_name); NTNFricRegNoRegularisation::dumpRestart(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegRubinAmpuero::readRestart(const std::string & file_name) { AKANTU_DEBUG_IN(); this->t_star.readRestartFile(file_name); NTNFricRegNoRegularisation::readRestart(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegRubinAmpuero::printself(std::ostream & stream, int indent) const { AKANTU_DEBUG_IN(); std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "NTNFricRegRubinAmpuero [" << std::endl; NTNFricRegNoRegularisation::printself(stream, ++indent); stream << space << "]" << std::endl; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegRubinAmpuero::addDumpFieldToDumper( const std::string & dumper_name, const std::string & field_id) { AKANTU_DEBUG_IN(); #ifdef AKANTU_USE_IOHELPER // const SynchronizedArray * nodal_filter = - // &(this->contact->getSlaves()); + // &(this->contact.getSlaves()); if (field_id == "t_star") { this->internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->t_star.getArray())); } else { NTNFricRegNoRegularisation::addDumpFieldToDumper(dumper_name, field_id); } #endif AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_rubin_ampuero.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_rubin_ampuero.hh index 371473983..953ac0a96 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_rubin_ampuero.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_rubin_ampuero.hh @@ -1,102 +1,102 @@ /** * @file ntn_fricreg_rubin_ampuero.hh * * @author David Simon Kammer * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Feb 23 2018 * * @brief regularisation that regularizes the contact pressure * * @section LICENSE * * Copyright (©) 2015-2018 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 __AST_NTN_FRICREG_RUBIN_AMPUERO_HH__ #define __AST_NTN_FRICREG_RUBIN_AMPUERO_HH__ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_fricreg_no_regularisation.hh" namespace akantu { /* -------------------------------------------------------------------------- */ class NTNFricRegRubinAmpuero : public NTNFricRegNoRegularisation { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - NTNFricRegRubinAmpuero(NTNBaseContact * contact, - const FrictionID & id = "rubin_ampuero", + NTNFricRegRubinAmpuero(NTNBaseContact & contact, + const ID & id = "rubin_ampuero", const MemoryID & memory_id = 0); virtual ~NTNFricRegRubinAmpuero(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual void registerSynchronizedArray(SynchronizedArrayBase & array); virtual void dumpRestart(const std::string & file_name) const; virtual void readRestart(const std::string & file_name); virtual void setToSteadyState(); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* Dumpable */ /* ------------------------------------------------------------------------ */ public: virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: protected: /// get the contact pressure (the norm: scalar value) virtual const SynchronizedArray & internalGetContactPressure(); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: SynchronizedArray t_star; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "ntn_fricreg_rubin_ampuero_inline_impl.cc" /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const NTNFricRegRubinAmpuero & _this) { _this.printself(stream); return stream; } } // namespace akantu #endif /* __AST_NTN_FRICREG_RUBIN_AMPUERO_HH__ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_simplified_prakash_clifton.cc b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_simplified_prakash_clifton.cc index 0c4bff531..74a8e7cf5 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_simplified_prakash_clifton.cc +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_simplified_prakash_clifton.cc @@ -1,163 +1,163 @@ /** * @file ntn_fricreg_simplified_prakash_clifton.cc * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief implementation of simplified prakash clifton with one parameter * * @section LICENSE * * Copyright (©) 2015-2018 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 . * */ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_fricreg_simplified_prakash_clifton.hh" #include "dumper_nodal_field.hh" #include "dumper_text.hh" namespace akantu { /* -------------------------------------------------------------------------- */ NTNFricRegSimplifiedPrakashClifton::NTNFricRegSimplifiedPrakashClifton( - NTNBaseContact * contact, const FrictionID & id, const MemoryID & memory_id) + NTNBaseContact & contact, const ID & id, const MemoryID & memory_id) : NTNFricRegNoRegularisation(contact, id, memory_id), t_star(0, 1, 0., id + ":t_star", 0., "t_star"), spc_internal(0, 1, 0., id + ":spc_internal", 0., "spc_internal") { AKANTU_DEBUG_IN(); NTNFricRegNoRegularisation::registerSynchronizedArray(this->t_star); NTNFricRegNoRegularisation::registerSynchronizedArray(this->spc_internal); this->registerParam("t_star", this->t_star, _pat_parsmod, "time scale of regularisation"); this->registerParam("spc_internal", this->spc_internal, _pat_internal, ""); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegSimplifiedPrakashClifton::computeFrictionalStrength() { AKANTU_DEBUG_IN(); - SolidMechanicsModel & model = this->contact->getModel(); + SolidMechanicsModel & model = this->contact.getModel(); Real delta_t = model.getTimeStep(); - UInt nb_contact_nodes = this->contact->getNbContactNodes(); + UInt nb_contact_nodes = this->contact.getNbContactNodes(); for (UInt n = 0; n < nb_contact_nodes; ++n) { Real alpha = delta_t / this->t_star(n); this->frictional_strength(n) += alpha * this->spc_internal(n); this->frictional_strength(n) /= 1 + alpha; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegSimplifiedPrakashClifton::setToSteadyState() { AKANTU_DEBUG_IN(); /// fill the spc_internal array computeFrictionalStrength(); /// set strength without regularisation - UInt nb_contact_nodes = this->contact->getNbContactNodes(); + UInt nb_contact_nodes = this->contact.getNbContactNodes(); for (UInt n = 0; n < nb_contact_nodes; ++n) { this->frictional_strength(n) = this->spc_internal(n); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegSimplifiedPrakashClifton::registerSynchronizedArray( SynchronizedArrayBase & array) { AKANTU_DEBUG_IN(); this->t_star.registerDependingArray(array); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegSimplifiedPrakashClifton::dumpRestart( const std::string & file_name) const { AKANTU_DEBUG_IN(); this->t_star.dumpRestartFile(file_name); this->spc_internal.dumpRestartFile(file_name); NTNFricRegNoRegularisation::dumpRestart(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegSimplifiedPrakashClifton::readRestart( const std::string & file_name) { AKANTU_DEBUG_IN(); this->t_star.readRestartFile(file_name); this->spc_internal.readRestartFile(file_name); NTNFricRegNoRegularisation::readRestart(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegSimplifiedPrakashClifton::printself(std::ostream & stream, int indent) const { AKANTU_DEBUG_IN(); std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "NTNFricRegSimplifiedPrakashClifton [" << std::endl; NTNFricRegNoRegularisation::printself(stream, ++indent); stream << space << "]" << std::endl; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNFricRegSimplifiedPrakashClifton::addDumpFieldToDumper( const std::string & dumper_name, const std::string & field_id) { AKANTU_DEBUG_IN(); #ifdef AKANTU_USE_IOHELPER // const SynchronizedArray * nodal_filter = - // &(this->contact->getSlaves()); + // &(this->contact.getSlaves()); if (field_id == "t_star") { this->internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->t_star.getArray())); } else { NTNFricRegNoRegularisation::addDumpFieldToDumper(dumper_name, field_id); } #endif AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_simplified_prakash_clifton.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_simplified_prakash_clifton.hh index b84fd0fe3..2d7519352 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_simplified_prakash_clifton.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/friction_regularisations/ntn_fricreg_simplified_prakash_clifton.hh @@ -1,114 +1,114 @@ /** * @file ntn_fricreg_simplified_prakash_clifton.hh * * @author David Simon Kammer * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Feb 23 2018 * * @brief regularisation that regularizes the frictional strength with one * parameter * * @section LICENSE * * Copyright (©) 2015-2018 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 __AST_NTN_FRICREG_SIMPLIFIED_PRAKASH_CLIFTON_HH__ #define __AST_NTN_FRICREG_SIMPLIFIED_PRAKASH_CLIFTON_HH__ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_fricreg_no_regularisation.hh" namespace akantu { /* -------------------------------------------------------------------------- */ class NTNFricRegSimplifiedPrakashClifton : public NTNFricRegNoRegularisation { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: NTNFricRegSimplifiedPrakashClifton( - NTNBaseContact * contact, - const FrictionID & id = "simplified_prakash_clifton", + NTNBaseContact & contact, + const ID & id = "simplified_prakash_clifton", const MemoryID & memory_id = 0); virtual ~NTNFricRegSimplifiedPrakashClifton(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual void registerSynchronizedArray(SynchronizedArrayBase & array); virtual void dumpRestart(const std::string & file_name) const; virtual void readRestart(const std::string & file_name); virtual void setToSteadyState(); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; protected: /// compute frictional strength according to friction law virtual void computeFrictionalStrength(); /* ------------------------------------------------------------------------ */ /* Dumpable */ /* ------------------------------------------------------------------------ */ public: virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: protected: /// get the frictional strength array virtual SynchronizedArray & internalGetFrictionalStrength() { return this->spc_internal; }; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: SynchronizedArray t_star; // to get the incremental frictional strength SynchronizedArray spc_internal; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "ntn_fricreg_simplified_prakash_clifton_inline_impl.cc" /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const NTNFricRegSimplifiedPrakashClifton & _this) { _this.printself(stream); return stream; } } // namespace akantu #endif /* __AST_NTN_FRICREG_SIMPLIFIED_PRAKASH_CLIFTON_HH__ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/mIIasym_contact.cc b/extra_packages/traction-at-split-node-contact/src/ntn_contact/mIIasym_contact.cc index e3608ce3a..76b19cdd8 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/mIIasym_contact.cc +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/mIIasym_contact.cc @@ -1,120 +1,120 @@ /** * @file mIIasym_contact.cc * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief * * @section LICENSE * * Copyright (©) 2015-2018 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 . * */ /* -------------------------------------------------------------------------- */ // simtools #include "mIIasym_contact.hh" namespace akantu { /* -------------------------------------------------------------------------- */ MIIASYMContact::MIIASYMContact(SolidMechanicsModel & model, - const ContactID & id, const MemoryID & memory_id) + const ID & id, const MemoryID & memory_id) : NTRFContact(model, id, memory_id) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MIIASYMContact::updateImpedance() { AKANTU_DEBUG_IN(); NTRFContact::updateImpedance(); for (UInt i = 0; i < this->impedance.size(); ++i) { this->impedance(i) *= 0.5; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /// WARNING: this is only valid for the acceleration in equilibrium void MIIASYMContact::computeRelativeNormalField( const Array & field, Array & rel_normal_field) const { AKANTU_DEBUG_IN(); NTRFContact::computeRelativeNormalField(field, rel_normal_field); for (auto it_rtfield = rel_normal_field.begin(); it_rtfield != rel_normal_field.end(); ++it_rtfield) { // in the anti-symmetric case *it_rtfield *= 2.; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MIIASYMContact::computeRelativeTangentialField( const Array & field, Array & rel_tang_field) const { AKANTU_DEBUG_IN(); NTRFContact::computeRelativeTangentialField(field, rel_tang_field); UInt dim = this->model.getSpatialDimension(); for (Array::iterator> it_rtfield = rel_tang_field.begin(dim); it_rtfield != rel_tang_field.end(dim); ++it_rtfield) { // in the anti-symmetric case, the tangential fields become twice as large *it_rtfield *= 2.; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MIIASYMContact::computeContactPressureInEquilibrium() { AKANTU_DEBUG_IN(); NTRFContact::computeContactPressure(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MIIASYMContact::printself(std::ostream & stream, int indent) const { AKANTU_DEBUG_IN(); std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "MIIASYMContact [" << std::endl; NTRFContact::printself(stream, indent); stream << space << "]" << std::endl; AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/mIIasym_contact.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/mIIasym_contact.hh index b2f8685ff..7e8d96be3 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/mIIasym_contact.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/mIIasym_contact.hh @@ -1,91 +1,91 @@ /** * @file mIIasym_contact.hh * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief contact for mode II anti-symmetric simulations * * @section LICENSE * * Copyright (©) 2015-2018 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 __AST_MIIASYM_CONTACT_HH__ #define __AST_MIIASYM_CONTACT_HH__ /* -------------------------------------------------------------------------- */ // simtools #include "ntrf_contact.hh" namespace akantu { /* -------------------------------------------------------------------------- */ class MIIASYMContact : public NTRFContact { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - MIIASYMContact(SolidMechanicsModel & model, const ContactID & id = "contact", + MIIASYMContact(SolidMechanicsModel & model, const ID & id = "contact", const MemoryID & memory_id = 0); - virtual ~MIIASYMContact(){}; + virtual ~MIIASYMContact() = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// update the impedance matrix virtual void updateImpedance(); /// compute contact pressure -> do nothing because can only compute it in /// equilibrium virtual void computeContactPressure(){}; /// compute relative normal field (only value that has to be multiplied with /// the normal) /// WARNING: this is only valid for the acceleration in equilibrium virtual void computeRelativeNormalField(const Array & field, Array & rel_normal_field) const; /// compute relative tangential field (complet array) /// relative to master nodes virtual void computeRelativeTangentialField(const Array & field, Array & rel_tang_field) const; /// compute contact pressure that is used over the entire time virtual void computeContactPressureInEquilibrium(); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const MIIASYMContact & _this) { _this.printself(stream); return stream; } } // namespace akantu #endif /* __AST_MIIASYM_CONTACT_HH__ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_contact.cc b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_contact.cc index 20b1eb424..f90facf28 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_contact.cc +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_contact.cc @@ -1,551 +1,555 @@ /** * @file ntn_base_contact.cc * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief implementation of ntn base contact * * @section LICENSE * * Copyright (©) 2015-2018 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 "ntn_base_contact.hh" +#include "dof_manager_default.hh" #include "dumpable_inline_impl.hh" #include "dumper_nodal_field.hh" #include "dumper_text.hh" #include "element_synchronizer.hh" #include "mesh_utils.hh" +#include "non_linear_solver_lumped.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ // NTNContactSynchElementFilter::NTNContactSynchElementFilter( -// NTNBaseContact * contact) +// NTNBaseContact & contact) // : contact(contact), -// connectivity(contact->getModel().getMesh().getConnectivities()) { +// connectivity(contact.getModel().getMesh().getConnectivities()) { // AKANTU_DEBUG_IN(); // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ // bool NTNContactSynchElementFilter::operator()(const Element & e) { // AKANTU_DEBUG_IN(); // ElementType type = e.type; // UInt element = e.element; // GhostType ghost_type = e.ghost_type; // // loop over all nodes of this element // bool need_element = false; // UInt nb_nodes = Mesh::getNbNodesPerElement(type); // for (UInt n = 0; n < nb_nodes; ++n) { // UInt nn = this->connectivity(type, ghost_type)(element, n); // // if one nodes is in this contact, we need this element -// if (this->contact->getNodeIndex(nn) >= 0) { +// if (this->contact.getNodeIndex(nn) >= 0) { // need_element = true; // break; // } // } // AKANTU_DEBUG_OUT(); // return need_element; // } /* -------------------------------------------------------------------------- */ NTNBaseContact::NTNBaseContact(SolidMechanicsModel & model, const ID & id, const MemoryID & memory_id) : Memory(id, memory_id), Dumpable(), model(model), slaves(0, 1, 0, id + ":slaves", std::numeric_limits::quiet_NaN(), "slaves"), normals(0, model.getSpatialDimension(), 0, id + ":normals", std::numeric_limits::quiet_NaN(), "normals"), contact_pressure( 0, model.getSpatialDimension(), 0, id + ":contact_pressure", std::numeric_limits::quiet_NaN(), "contact_pressure"), is_in_contact(0, 1, false, id + ":is_in_contact", false, "is_in_contact"), lumped_boundary_slaves(0, 1, 0, id + ":lumped_boundary_slaves", std::numeric_limits::quiet_NaN(), "lumped_boundary_slaves"), impedance(0, 1, 0, id + ":impedance", std::numeric_limits::quiet_NaN(), "impedance"), contact_surfaces(), slave_elements("slave_elements", id, memory_id), node_to_elements() { AKANTU_DEBUG_IN(); FEEngine & boundary_fem = this->model.getFEEngineBoundary(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { boundary_fem.initShapeFunctions(*gt); } Mesh & mesh = this->model.getMesh(); UInt spatial_dimension = this->model.getSpatialDimension(); this->slave_elements.initialize(mesh, _spatial_dimension = spatial_dimension - 1); MeshUtils::buildNode2Elements(mesh, this->node_to_elements, spatial_dimension - 1); this->registerDumper("text_all", id, true); this->addDumpFilteredMesh(mesh, slave_elements, slaves.getArray(), spatial_dimension - 1, _not_ghost, _ek_regular); // parallelisation this->synch_registry = std::make_unique(); this->synch_registry->registerDataAccessor(*this); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ NTNBaseContact::~NTNBaseContact() = default; /* -------------------------------------------------------------------------- */ void NTNBaseContact::initParallel() { AKANTU_DEBUG_IN(); this->synchronizer = std::make_unique( this->model.getMesh().getElementSynchronizer()); this->synchronizer->filterScheme([&](auto && element) { // loop over all nodes of this element Vector conn = const_cast(this->model.getMesh()) .getConnectivity(element); for (auto & node : conn) { // if one nodes is in this contact, we need this element if (this->getNodeIndex(node) >= 0) { return true; } } return false; }); this->synch_registry->registerSynchronizer(*(this->synchronizer), _gst_cf_nodal); this->synch_registry->registerSynchronizer(*(this->synchronizer), _gst_cf_incr); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::findBoundaryElements( const Array & interface_nodes, ElementTypeMapArray & elements) { AKANTU_DEBUG_IN(); // add connected boundary elements that have all nodes on this contact for (const auto & node : interface_nodes) { for (const auto & element : this->node_to_elements.getRow(node)) { Vector conn = const_cast(this->model.getMesh()) .getConnectivity(element); auto nb_nodes = conn.size(); decltype(nb_nodes) nb_found_nodes = 0; for (auto & nn : conn) { if (interface_nodes.find(nn) != UInt(-1)) { nb_found_nodes++; } else { break; } } // this is an element between all contact nodes // and is not already in the elements if ((nb_found_nodes == nb_nodes) && (elements(element.type, element.ghost_type).find(element.element) == UInt(-1))) { elements(element.type, element.ghost_type).push_back(element.element); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::addSplitNode(UInt node) { AKANTU_DEBUG_IN(); UInt dim = this->model.getSpatialDimension(); // add to node arrays this->slaves.push_back(node); // set contact as false this->is_in_contact.push_back(false); // before initializing // set contact pressure, normal, lumped_boundary to Nan this->contact_pressure.push_back(std::numeric_limits::quiet_NaN()); this->impedance.push_back(std::numeric_limits::quiet_NaN()); this->lumped_boundary_slaves.push_back( std::numeric_limits::quiet_NaN()); Vector nan_normal(dim, std::numeric_limits::quiet_NaN()); this->normals.push_back(nan_normal); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::registerSynchronizedArray(SynchronizedArrayBase & array) { AKANTU_DEBUG_IN(); this->slaves.registerDependingArray(array); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::dumpRestart(const std::string & file_name) const { AKANTU_DEBUG_IN(); this->slaves.dumpRestartFile(file_name); this->normals.dumpRestartFile(file_name); this->is_in_contact.dumpRestartFile(file_name); this->contact_pressure.dumpRestartFile(file_name); this->lumped_boundary_slaves.dumpRestartFile(file_name); this->impedance.dumpRestartFile(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::readRestart(const std::string & file_name) { AKANTU_DEBUG_IN(); this->slaves.readRestartFile(file_name); this->normals.readRestartFile(file_name); this->is_in_contact.readRestartFile(file_name); this->contact_pressure.readRestartFile(file_name); this->lumped_boundary_slaves.readRestartFile(file_name); this->impedance.readRestartFile(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ UInt NTNBaseContact::getNbNodesInContact() const { AKANTU_DEBUG_IN(); UInt nb_contact = 0; UInt nb_nodes = this->getNbContactNodes(); const Mesh & mesh = this->model.getMesh(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(this->slaves(n)); - bool is_pbc_slave_node = this->model.isPBCSlaveNode(this->slaves(n)); + bool is_pbc_slave_node = mesh.isPeriodicSlave(this->slaves(n)); if (is_local_node && !is_pbc_slave_node && this->is_in_contact(n)) { nb_contact++; } } mesh.getCommunicator().allReduce(nb_contact, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return nb_contact; } /* -------------------------------------------------------------------------- */ void NTNBaseContact::updateInternalData() { AKANTU_DEBUG_IN(); updateNormals(); updateLumpedBoundary(); updateImpedance(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::updateLumpedBoundary() { AKANTU_DEBUG_IN(); this->internalUpdateLumpedBoundary(this->slaves.getArray(), this->slave_elements, this->lumped_boundary_slaves); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::internalUpdateLumpedBoundary( const Array & nodes, const ElementTypeMapArray & elements, SynchronizedArray & boundary) { AKANTU_DEBUG_IN(); // set all values in lumped_boundary to zero boundary.clear(); UInt dim = this->model.getSpatialDimension(); // UInt nb_contact_nodes = getNbContactNodes(); const FEEngine & boundary_fem = this->model.getFEEngineBoundary(); const Mesh & mesh = this->model.getMesh(); - for (ghost_type_t::iterator gt = ghost_type_t::begin(); - gt != ghost_type_t::end(); ++gt) { - Mesh::type_iterator it = mesh.firstType(dim - 1, *gt); - Mesh::type_iterator last = mesh.lastType(dim - 1, *gt); - for (; it != last; ++it) { - UInt nb_elements = mesh.getNbElement(*it, *gt); - UInt nb_nodes_per_element = mesh.getNbNodesPerElement(*it); - const Array & connectivity = mesh.getConnectivity(*it, *gt); + for(auto ghost_type : ghost_types) { + for (auto & type : mesh.elementTypes(dim - 1, ghost_type)) { + UInt nb_elements = mesh.getNbElement(type, ghost_type); + UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); + const Array & connectivity = mesh.getConnectivity(type, ghost_type); // get shapes and compute integral - const Array & shapes = boundary_fem.getShapes(*it, *gt); + const Array & shapes = boundary_fem.getShapes(type, ghost_type); Array area(nb_elements, nb_nodes_per_element); - boundary_fem.integrate(shapes, area, nb_nodes_per_element, *it, *gt); + boundary_fem.integrate(shapes, area, nb_nodes_per_element, type, ghost_type); if (this->contact_surfaces.size() == 0) { AKANTU_DEBUG_WARNING( "No surfaces in ntn base contact." << " You have to define the lumped boundary by yourself."); } - Array::const_iterator elem_it = (elements)(*it, *gt).begin(); + Array::const_iterator elem_it = (elements)(type, ghost_type).begin(); Array::const_iterator elem_it_end = - (elements)(*it, *gt).end(); + (elements)(type, ghost_type).end(); // loop over contact nodes for (; elem_it != elem_it_end; ++elem_it) { for (UInt q = 0; q < nb_nodes_per_element; ++q) { UInt node = connectivity(*elem_it, q); UInt node_index = nodes.find(node); - AKANTU_DEBUG_ASSERT(node_index != UInt(-1), - "Could not find node " << node - << " in the array!"); + AKANTU_DEBUG_ASSERT(node_index != UInt(-1), "Could not find node " + << node + << " in the array!"); Real area_to_add = area(*elem_it, q); boundary(node_index) += area_to_add; } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::computeContactPressure() { AKANTU_DEBUG_IN(); UInt dim = this->model.getSpatialDimension(); Real delta_t = this->model.getTimeStep(); UInt nb_contact_nodes = getNbContactNodes(); AKANTU_DEBUG_ASSERT(delta_t > 0., "Cannot compute contact pressure if no time step is set"); // synchronize data this->synch_registry->synchronize(_gst_cf_nodal); + auto && dof_manager = + dynamic_cast(model.getDOFManager()); + const auto & b = dof_manager.getResidual(); + Array acceleration(b.size(), dim); + const auto & blocked_dofs = dof_manager.getGlobalBlockedDOFs(); + const auto & A = dof_manager.getLumpedMatrix("M"); + // pre-compute the acceleration // (not increment acceleration, because residual is still Kf) - Array acceleration(this->model.getMesh().getNbNodes(), dim, 0.); - this->model.solveLumped(acceleration, this->model.getMass(), - this->model.getInternalForce(), - this->model.getBlockedDOFs(), this->model.getF_M2A()); + NonLinearSolverLumped::solveLumped(A, acceleration, b, blocked_dofs, + this->model.getF_M2A()); // compute relative normal fields of displacement, velocity and acceleration Array r_disp(0, 1); Array r_velo(0, 1); Array r_acce(0, 1); Array r_old_acce(0, 1); computeNormalGap(r_disp); // computeRelativeNormalField(this->model.getCurrentPosition(), r_disp); computeRelativeNormalField(this->model.getVelocity(), r_velo); computeRelativeNormalField(acceleration, r_acce); computeRelativeNormalField(this->model.getAcceleration(), r_old_acce); AKANTU_DEBUG_ASSERT(r_disp.size() == nb_contact_nodes, "computeRelativeNormalField does not give back arrays " << "size == nb_contact_nodes. nb_contact_nodes = " << nb_contact_nodes << " | array size = " << r_disp.size()); // compute gap array for all nodes Array gap(nb_contact_nodes, 1); Real * gap_p = gap.storage(); Real * r_disp_p = r_disp.storage(); Real * r_velo_p = r_velo.storage(); Real * r_acce_p = r_acce.storage(); Real * r_old_acce_p = r_old_acce.storage(); for (UInt i = 0; i < nb_contact_nodes; ++i) { *gap_p = *r_disp_p + delta_t * *r_velo_p + delta_t * delta_t * *r_acce_p - 0.5 * delta_t * delta_t * *r_old_acce_p; // increment pointers gap_p++; r_disp_p++; r_velo_p++; r_acce_p++; r_old_acce_p++; } // check if gap is negative -> is in contact for (UInt n = 0; n < nb_contact_nodes; ++n) { if (gap(n) <= 0.) { for (UInt d = 0; d < dim; ++d) { this->contact_pressure(n, d) = this->impedance(n) * gap(n) / (2 * delta_t) * this->normals(n, d); } this->is_in_contact(n) = true; } else { for (UInt d = 0; d < dim; ++d) this->contact_pressure(n, d) = 0.; this->is_in_contact(n) = false; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::applyContactPressure() { AKANTU_DEBUG_IN(); UInt nb_contact_nodes = getNbContactNodes(); UInt dim = this->model.getSpatialDimension(); Array & residual = this->model.getInternalForce(); for (UInt n = 0; n < nb_contact_nodes; ++n) { UInt slave = this->slaves(n); for (UInt d = 0; d < dim; ++d) { // residual(master,d) += this->lumped_boundary(n,0) * // this->contact_pressure(n,d); residual(slave, d) -= this->lumped_boundary_slaves(n) * this->contact_pressure(n, d); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Int NTNBaseContact::getNodeIndex(UInt node) const { return this->slaves.find(node); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::printself(std::ostream & stream, int indent) const { AKANTU_DEBUG_IN(); std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "NTNBaseContact [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + slaves : " << std::endl; this->slaves.printself(stream, indent + 2); stream << space << " + normals : " << std::endl; this->normals.printself(stream, indent + 2); stream << space << " + contact_pressure : " << std::endl; this->contact_pressure.printself(stream, indent + 2); stream << space << "]" << std::endl; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::syncArrays(SyncChoice sync_choice) { AKANTU_DEBUG_IN(); this->slaves.syncElements(sync_choice); this->normals.syncElements(sync_choice); this->is_in_contact.syncElements(sync_choice); this->contact_pressure.syncElements(sync_choice); this->lumped_boundary_slaves.syncElements(sync_choice); this->impedance.syncElements(sync_choice); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseContact::addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id) { AKANTU_DEBUG_IN(); #ifdef AKANTU_USE_IOHELPER const Array & nodal_filter = this->slaves.getArray(); #define ADD_FIELD(field_id, field, type) \ internalAddDumpFieldToDumper( \ dumper_name, field_id, \ new dumper::NodalField, Array>( \ field, 0, 0, &nodal_filter)) if (field_id == "displacement") { ADD_FIELD(field_id, this->model.getDisplacement(), Real); } else if (field_id == "mass") { ADD_FIELD(field_id, this->model.getMass(), Real); } else if (field_id == "velocity") { ADD_FIELD(field_id, this->model.getVelocity(), Real); } else if (field_id == "acceleration") { ADD_FIELD(field_id, this->model.getAcceleration(), Real); } else if (field_id == "external_force") { ADD_FIELD(field_id, this->model.getForce(), Real); } else if (field_id == "internal_force") { ADD_FIELD(field_id, this->model.getInternalForce(), Real); } else if (field_id == "blocked_dofs") { ADD_FIELD(field_id, this->model.getBlockedDOFs(), bool); } else if (field_id == "increment") { ADD_FIELD(field_id, this->model.getIncrement(), Real); } else if (field_id == "normal") { internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->normals.getArray())); } else if (field_id == "contact_pressure") { internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->contact_pressure.getArray())); } else if (field_id == "is_in_contact") { internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->is_in_contact.getArray())); } else if (field_id == "lumped_boundary_slave") { internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->lumped_boundary_slaves.getArray())); } else if (field_id == "impedance") { internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->impedance.getArray())); } else { std::cerr << "Could not add field '" << field_id << "' to the dumper. Just ignored it." << std::endl; } #undef ADD_FIELD #endif AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_contact.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_contact.hh index 407ab5f73..347375c35 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_contact.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_contact.hh @@ -1,246 +1,250 @@ /** * @file ntn_base_contact.hh * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief base contact for ntn and ntrf contact * * @section LICENSE * * Copyright (©) 2015-2018 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 __AST_NTN_BASE_CONTACT_HH__ #define __AST_NTN_BASE_CONTACT_HH__ /* -------------------------------------------------------------------------- */ // akantu #include "aka_csr.hh" #include "solid_mechanics_model.hh" // simtools #include "synchronized_array.hh" namespace akantu { class NTNBaseContact; /* -------------------------------------------------------------------------- */ // class NTNContactSynchElementFilter : public SynchElementFilter { // public: // // constructor -// NTNContactSynchElementFilter(NTNBaseContact * contact); +// NTNContactSynchElementFilter(NTNBaseContact & contact); // // answer to: do we need this element ? // virtual bool operator()(const Element & e); // private: -// const NTNBaseContact * contact; +// const NTNBaseContact & contact; // const ElementTypeMapArray & connectivity; // }; /* -------------------------------------------------------------------------- */ class NTNBaseContact : protected Memory, public DataAccessor, public Dumpable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: NTNBaseContact(SolidMechanicsModel & model, const ID & id = "contact", const MemoryID & memory_id = 0); virtual ~NTNBaseContact(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initializes ntn contact parallel virtual void initParallel(); /// add split node virtual void addSplitNode(UInt node); /// update normals, lumped boundary, and impedance virtual void updateInternalData(); /// update (compute the normals) virtual void updateNormals() = 0; /// update the lumped boundary B matrix virtual void updateLumpedBoundary(); /// update the impedance matrix virtual void updateImpedance() = 0; /// compute the normal contact force virtual void computeContactPressure(); /// impose the normal contact force virtual void applyContactPressure(); /// register synchronizedarrays for sync virtual void registerSynchronizedArray(SynchronizedArrayBase & array); /// dump restart file virtual void dumpRestart(const std::string & file_name) const; /// read restart file virtual void readRestart(const std::string & file_name); /// compute the normal gap virtual void computeNormalGap(Array & gap) const = 0; /// compute relative normal field (only value that has to be multiplied with /// the normal) /// relative to master nodes virtual void computeRelativeNormalField(const Array & field, Array & rel_normal_field) const = 0; /// compute relative tangential field (complet array) /// relative to master nodes virtual void computeRelativeTangentialField(const Array & field, Array & rel_tang_field) const = 0; /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; protected: /// updateLumpedBoundary virtual void internalUpdateLumpedBoundary(const Array & nodes, const ElementTypeMapArray & elements, SynchronizedArray & boundary); // to find the slave_elements or master_elements virtual void findBoundaryElements(const Array & interface_nodes, ElementTypeMapArray & elements); /// synchronize arrays virtual void syncArrays(SyncChoice sync_choice); /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: inline UInt getNbData(const Array & elements, const SynchronizationTag & tag) const override; inline void packData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) const override; inline void unpackData(CommunicationBuffer & buffer, const Array & elements, const SynchronizationTag & tag) override; /* ------------------------------------------------------------------------ */ /* Dumpable */ /* ------------------------------------------------------------------------ */ public: virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Model, model, SolidMechanicsModel &) AKANTU_GET_MACRO(Slaves, slaves, const SynchronizedArray &) AKANTU_GET_MACRO(Normals, normals, const SynchronizedArray &) AKANTU_GET_MACRO(ContactPressure, contact_pressure, const SynchronizedArray &) AKANTU_GET_MACRO(LumpedBoundarySlaves, lumped_boundary_slaves, const SynchronizedArray &) AKANTU_GET_MACRO(Impedance, impedance, const SynchronizedArray &) AKANTU_GET_MACRO(IsInContact, is_in_contact, const SynchronizedArray &) AKANTU_GET_MACRO(SlaveElements, slave_elements, const ElementTypeMapArray &) AKANTU_GET_MACRO(SynchronizerRegistry, *synch_registry, SynchronizerRegistry &) /// get number of nodes that are in contact (globally, on all procs together) /// is_in_contact = true virtual UInt getNbNodesInContact() const; /// get index of node in either slaves or masters array /// if node is in neither of them, return -1 virtual Int getNodeIndex(UInt node) const; /// get number of contact nodes: nodes in the system locally (on this proc) /// is_in_contact = true and false, because just in the system - virtual UInt getNbContactNodes() const { return this->slaves.size(); }; + virtual UInt getNbContactNodes() const { return this->slaves.size(); } + + bool isNTNContact() const { return this->is_ntn_contact; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: typedef std::set SurfacePtrSet; SolidMechanicsModel & model; /// array of slave nodes SynchronizedArray slaves; /// array of normals SynchronizedArray normals; /// array indicating if nodes are in contact SynchronizedArray contact_pressure; /// array indicating if nodes are in contact SynchronizedArray is_in_contact; /// boundary matrix for slave nodes SynchronizedArray lumped_boundary_slaves; /// impedance matrix SynchronizedArray impedance; /// contact surface SurfacePtrSet contact_surfaces; /// element list for dump and lumped_boundary ElementTypeMapArray slave_elements; CSR node_to_elements; /// parallelisation std::unique_ptr synch_registry; std::unique_ptr synchronizer; + + bool is_ntn_contact{true}; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "ntn_base_contact_inline_impl.cc" /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const NTNBaseContact & _this) { _this.printself(stream); return stream; } } // namespace akantu #endif /* __AST_NTN_BASE_CONTACT_HH__ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_friction.cc b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_friction.cc index df984c08d..fe39d29be 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_friction.cc +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_friction.cc @@ -1,382 +1,383 @@ /** * @file ntn_base_friction.cc * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief implementation of ntn base friction * * @section LICENSE * * Copyright (©) 2015-2018 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 . * */ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_base_friction.hh" +#include "dof_manager_default.hh" #include "dumper_nodal_field.hh" #include "dumper_text.hh" +#include "non_linear_solver_lumped.hh" namespace akantu { /* -------------------------------------------------------------------------- */ -NTNBaseFriction::NTNBaseFriction(NTNBaseContact * contact, - const FrictionID & id, +NTNBaseFriction::NTNBaseFriction(NTNBaseContact & contact, const ID & id, const MemoryID & memory_id) - : Memory(id, memory_id), Parsable(_st_friction, id), Dumpable(), + : Memory(id, memory_id), Parsable(ParserType::_friction, id), Dumpable(), contact(contact), is_sticking(0, 1, true, id + ":is_sticking", true, "is_sticking"), frictional_strength(0, 1, 0., id + ":frictional_strength", 0., "frictional_strength"), - friction_traction(0, contact->getModel().getSpatialDimension(), 0., + friction_traction(0, contact.getModel().getSpatialDimension(), 0., id + ":friction_traction", 0., "friction_traction"), slip(0, 1, 0., id + ":slip", 0., "slip"), cumulative_slip(0, 1, 0., id + ":cumulative_slip", 0., "cumulative_slip"), - slip_velocity(0, contact->getModel().getSpatialDimension(), 0., + slip_velocity(0, contact.getModel().getSpatialDimension(), 0., id + ":slip_velocity", 0., "slip_velocity") { AKANTU_DEBUG_IN(); - this->contact->registerSynchronizedArray(this->is_sticking); - this->contact->registerSynchronizedArray(this->frictional_strength); - this->contact->registerSynchronizedArray(this->friction_traction); - this->contact->registerSynchronizedArray(this->slip); - this->contact->registerSynchronizedArray(this->cumulative_slip); - this->contact->registerSynchronizedArray(this->slip_velocity); + this->contact.registerSynchronizedArray(this->is_sticking); + this->contact.registerSynchronizedArray(this->frictional_strength); + this->contact.registerSynchronizedArray(this->friction_traction); + this->contact.registerSynchronizedArray(this->slip); + this->contact.registerSynchronizedArray(this->cumulative_slip); + this->contact.registerSynchronizedArray(this->slip_velocity); - contact->getModel().setIncrementFlagOn(); - - this->registerExternalDumper(contact->getDumper(), - contact->getDefaultDumperName(), true); + this->registerExternalDumper(contact.getDumper(), + contact.getDefaultDumperName(), true); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseFriction::updateSlip() { AKANTU_DEBUG_IN(); - SolidMechanicsModel & model = this->contact->getModel(); + SolidMechanicsModel & model = this->contact.getModel(); UInt dim = model.getSpatialDimension(); // synchronize increment - this->contact->getSynchronizerRegistry()->synchronize(_gst_cf_incr); + this->contact.getSynchronizerRegistry().synchronize(_gst_cf_incr); Array rel_tan_incr(0, dim); - this->contact->computeRelativeTangentialField(model.getIncrement(), + this->contact.computeRelativeTangentialField(model.getIncrement(), rel_tan_incr); Array::const_iterator> it = rel_tan_incr.begin(dim); - UInt nb_nodes = this->contact->getNbContactNodes(); + UInt nb_nodes = this->contact.getNbContactNodes(); for (UInt n = 0; n < nb_nodes; ++n) { if (this->is_sticking(n)) { this->slip(n) = 0.; } else { const Vector & rti = it[n]; this->slip(n) += rti.norm(); this->cumulative_slip(n) += rti.norm(); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseFriction::computeFrictionTraction() { AKANTU_DEBUG_IN(); this->computeStickTraction(); this->computeFrictionalStrength(); - SolidMechanicsModel & model = this->contact->getModel(); + SolidMechanicsModel & model = this->contact.getModel(); UInt dim = model.getSpatialDimension(); // get contact arrays const SynchronizedArray & is_in_contact = - this->contact->getIsInContact(); + this->contact.getIsInContact(); Array & traction = const_cast &>(this->friction_traction.getArray()); Array::iterator> it_fric_trac = traction.begin(dim); this->is_sticking.clear(); // set to not sticking - UInt nb_contact_nodes = this->contact->getNbContactNodes(); + UInt nb_contact_nodes = this->contact.getNbContactNodes(); for (UInt n = 0; n < nb_contact_nodes; ++n) { // node pair is in contact if (is_in_contact(n)) { Vector fric_trac = it_fric_trac[n]; // check if it is larger than frictional strength Real abs_fric = fric_trac.norm(); if (abs_fric != 0.) { Real alpha = this->frictional_strength(n) / abs_fric; // larger -> sliding if (alpha < 1.) { fric_trac *= alpha; } else this->is_sticking(n) = true; } else { // frictional traction is already zero this->is_sticking(n) = true; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseFriction::computeStickTraction() { AKANTU_DEBUG_IN(); - SolidMechanicsModel & model = this->contact->getModel(); + SolidMechanicsModel & model = this->contact.getModel(); UInt dim = model.getSpatialDimension(); Real delta_t = model.getTimeStep(); - UInt nb_contact_nodes = this->contact->getNbContactNodes(); + UInt nb_contact_nodes = this->contact.getNbContactNodes(); // get contact arrays - const SynchronizedArray & impedance = this->contact->getImpedance(); + const SynchronizedArray & impedance = this->contact.getImpedance(); const SynchronizedArray & is_in_contact = - this->contact->getIsInContact(); + this->contact.getIsInContact(); + + auto && dof_manager = + dynamic_cast(model.getDOFManager()); + const auto & b = dof_manager.getResidual(); + Array acceleration(b.size(), dim); + const auto & blocked_dofs = dof_manager.getGlobalBlockedDOFs(); + const auto & A = dof_manager.getLumpedMatrix("M"); // pre-compute the acceleration // (not increment acceleration, because residual is still Kf) - Array acceleration(model.getFEEngine().getMesh().getNbNodes(), dim, 0.); - model.solveLumped(acceleration, model.getMass(), model.getResidual(), - model.getBlockedDOFs(), model.getF_M2A()); + NonLinearSolverLumped::solveLumped(A, acceleration, b, blocked_dofs, + model.getF_M2A()); // compute relative normal fields of velocity and acceleration Array r_velo(0, dim); Array r_acce(0, dim); Array r_old_acce(0, dim); - this->contact->computeRelativeTangentialField(model.getVelocity(), r_velo); - this->contact->computeRelativeTangentialField(acceleration, r_acce); - this->contact->computeRelativeTangentialField(model.getAcceleration(), + this->contact.computeRelativeTangentialField(model.getVelocity(), r_velo); + this->contact.computeRelativeTangentialField(acceleration, r_acce); + this->contact.computeRelativeTangentialField(model.getAcceleration(), r_old_acce); AKANTU_DEBUG_ASSERT(r_velo.size() == nb_contact_nodes, "computeRelativeNormalField does not give back arrays " << "size == nb_contact_nodes. nb_contact_nodes = " << nb_contact_nodes << " | array size = " << r_velo.size()); // compute tangential gap_dot array for all nodes Array gap_dot(nb_contact_nodes, dim); - Real * gap_dot_p = gap_dot.storage(); - Real * r_velo_p = r_velo.storage(); - Real * r_acce_p = r_acce.storage(); - Real * r_old_acce_p = r_old_acce.storage(); - for (UInt i = 0; i < nb_contact_nodes * dim; ++i) { - *gap_dot_p = - *r_velo_p + delta_t * *r_acce_p - 0.5 * delta_t * *r_old_acce_p; - // increment pointers - gap_dot_p++; - r_velo_p++; - r_acce_p++; - r_old_acce_p++; + for (auto && data : zip(make_view(gap_dot), make_view(r_velo), + make_view(r_acce), make_view(r_old_acce))) { + auto & gap_dot = std::get<0>(data); + auto & r_velo = std::get<1>(data); + auto & r_acce = std::get<2>(data); + auto & r_old_acce = std::get<3>(data); + + gap_dot = r_velo + delta_t * r_acce - 1. / 2. * delta_t * r_old_acce; } // compute friction traction to stop sliding Array & traction = const_cast &>(this->friction_traction.getArray()); auto it_fric_trac = traction.begin(dim); for (UInt n = 0; n < nb_contact_nodes; ++n) { Vector fric_trac = it_fric_trac[n]; // node pair is NOT in contact if (!is_in_contact(n)) { fric_trac.clear(); // set to zero } // node pair is in contact else { // compute friction traction for (UInt d = 0; d < dim; ++d) fric_trac(d) = impedance(n) * gap_dot(n, d) / 2.; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseFriction::applyFrictionTraction() { AKANTU_DEBUG_IN(); - SolidMechanicsModel & model = this->contact->getModel(); - Array & residual = model.getResidual(); + SolidMechanicsModel & model = this->contact.getModel(); + Array & residual = model.getInternalForce(); UInt dim = model.getSpatialDimension(); - const SynchronizedArray & slaves = this->contact->getSlaves(); + const SynchronizedArray & slaves = this->contact.getSlaves(); const SynchronizedArray & lumped_boundary_slaves = - this->contact->getLumpedBoundarySlaves(); + this->contact.getLumpedBoundarySlaves(); - UInt nb_contact_nodes = this->contact->getNbContactNodes(); + UInt nb_contact_nodes = this->contact.getNbContactNodes(); for (UInt n = 0; n < nb_contact_nodes; ++n) { UInt slave = slaves(n); for (UInt d = 0; d < dim; ++d) { residual(slave, d) -= lumped_boundary_slaves(n) * this->friction_traction(n, d); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseFriction::registerSynchronizedArray(SynchronizedArrayBase & array) { AKANTU_DEBUG_IN(); this->frictional_strength.registerDependingArray(array); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseFriction::dumpRestart(const std::string & file_name) const { AKANTU_DEBUG_IN(); this->is_sticking.dumpRestartFile(file_name); this->frictional_strength.dumpRestartFile(file_name); this->friction_traction.dumpRestartFile(file_name); this->slip.dumpRestartFile(file_name); this->cumulative_slip.dumpRestartFile(file_name); this->slip_velocity.dumpRestartFile(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseFriction::readRestart(const std::string & file_name) { AKANTU_DEBUG_IN(); this->is_sticking.readRestartFile(file_name); this->frictional_strength.readRestartFile(file_name); this->friction_traction.readRestartFile(file_name); this->cumulative_slip.readRestartFile(file_name); this->slip_velocity.readRestartFile(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseFriction::setParam(const std::string & name, UInt node, Real value) { AKANTU_DEBUG_IN(); - SynchronizedArray & array = this->get>(name); - Int index = this->contact->getNodeIndex(node); + SynchronizedArray & array = this->get(name).get>(); + Int index = this->contact.getNodeIndex(node); if (index < 0) { AKANTU_DEBUG_WARNING("Node " << node << " is not a contact node. " << "Therefore, cannot set interface parameter!!"); } else { array(index) = value; // put value } AKANTU_DEBUG_OUT(); -}; +} /* -------------------------------------------------------------------------- */ UInt NTNBaseFriction::getNbStickingNodes() const { AKANTU_DEBUG_IN(); UInt nb_stick = 0; - UInt nb_nodes = this->contact->getNbContactNodes(); - const SynchronizedArray & nodes = this->contact->getSlaves(); + UInt nb_nodes = this->contact.getNbContactNodes(); + const SynchronizedArray & nodes = this->contact.getSlaves(); const SynchronizedArray & is_in_contact = - this->contact->getIsInContact(); + this->contact.getIsInContact(); - const Mesh & mesh = this->contact->getModel().getMesh(); + const Mesh & mesh = this->contact.getModel().getMesh(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(nodes(n)); - bool is_pbc_slave_node = this->contact->getModel().isPBCSlaveNode(nodes(n)); + bool is_pbc_slave_node = mesh.isPeriodicSlave(nodes(n)); if (is_local_node && !is_pbc_slave_node && is_in_contact(n) && this->is_sticking(n)) { nb_stick++; } } - StaticCommunicator::getStaticCommunicator().allReduce(&nb_stick, 1, _so_sum); + mesh.getCommunicator().allReduce(nb_stick, SynchronizerOperation::_sum); AKANTU_DEBUG_OUT(); return nb_stick; } /* -------------------------------------------------------------------------- */ void NTNBaseFriction::printself(std::ostream & stream, int indent) const { AKANTU_DEBUG_IN(); std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "NTNBaseFriction [" << std::endl; Parsable::printself(stream, indent); stream << space << "]" << std::endl; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNBaseFriction::addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id) { AKANTU_DEBUG_IN(); #ifdef AKANTU_USE_IOHELPER // const SynchronizedArray * nodal_filter = - // &(this->contact->getSlaves()); + // &(this->contact.getSlaves()); if (field_id == "is_sticking") { this->internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->is_sticking.getArray())); } else if (field_id == "frictional_strength") { this->internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->frictional_strength.getArray())); } else if (field_id == "friction_traction") { this->internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->friction_traction.getArray())); } else if (field_id == "slip") { this->internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->slip.getArray())); } else if (field_id == "cumulative_slip") { this->internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->cumulative_slip.getArray())); } else if (field_id == "slip_velocity") { this->internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->slip_velocity.getArray())); } else { - this->contact->addDumpFieldToDumper(dumper_name, field_id); + this->contact.addDumpFieldToDumper(dumper_name, field_id); } #endif AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_friction.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_friction.hh index a2e7ca7ae..6053b29eb 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_friction.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_base_friction.hh @@ -1,177 +1,179 @@ /** * @file ntn_base_friction.hh * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief base class for ntn and ntrf friction * * @section LICENSE * * Copyright (©) 2015-2018 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 __AST_NTN_BASE_FRICTION_HH__ #define __AST_NTN_BASE_FRICTION_HH__ /* -------------------------------------------------------------------------- */ // akantu #include "parsable.hh" // simtools #include "ntn_base_contact.hh" /* -------------------------------------------------------------------------- */ namespace akantu { + +/* -------------------------------------------------------------------------- */ template <> -inline void ParsableParamTyped>::parseParam( - const ParserParameter & in_param) { - ParsableParam::parseParam(in_param); - Real tmp = in_param; - param.setAndChangeDefault(tmp); +inline void +ParameterTyped>::setAuto(const ParserParameter & in_param) { + Parameter::setAuto(in_param); + Real r = in_param; + param.setAndChangeDefault(r); } /* -------------------------------------------------------------------------- */ template <> template <> -inline void ParsableParamTyped>::setTyped( +inline void ParameterTyped>::setTyped( const Real & value) { param.setAndChangeDefault(value); } /* -------------------------------------------------------------------------- */ class NTNBaseFriction : protected Memory, public Parsable, public Dumpable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - NTNBaseFriction(NTNBaseContact * contact, const FrictionID & id = "friction", + NTNBaseFriction(NTNBaseContact & contact, const ID & id = "friction", const MemoryID & memory_id = 0); - virtual ~NTNBaseFriction(){}; + virtual ~NTNBaseFriction() = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// compute friction traction virtual void computeFrictionTraction(); /// compute stick traction (friction traction needed to stick the nodes) virtual void computeStickTraction(); /// apply the friction force virtual void applyFrictionTraction(); /// compute slip virtual void updateSlip(); /// register Syncronizedarrays for sync virtual void registerSynchronizedArray(SynchronizedArrayBase & array); /// dump restart file virtual void dumpRestart(const std::string & file_name) const; /// read restart file virtual void readRestart(const std::string & file_name); /// set to steady state virtual void setToSteadyState() { AKANTU_TO_IMPLEMENT(); }; /// get the number of sticking nodes (in parallel) /// a node that is not in contact does not count as sticking virtual UInt getNbStickingNodes() const; /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; protected: /// compute frictional strength according to friction law virtual void computeFrictionalStrength() = 0; /* ------------------------------------------------------------------------ */ /* Dumpable */ /* ------------------------------------------------------------------------ */ public: virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: - AKANTU_GET_MACRO(Contact, contact, const NTNBaseContact *) + AKANTU_GET_MACRO(Contact, contact, const NTNBaseContact &) AKANTU_GET_MACRO(IsSticking, is_sticking, const SynchronizedArray &) AKANTU_GET_MACRO(FrictionalStrength, frictional_strength, const SynchronizedArray &) AKANTU_GET_MACRO(FrictionTraction, friction_traction, const SynchronizedArray &) AKANTU_GET_MACRO(Slip, slip, const SynchronizedArray &) AKANTU_GET_MACRO(CumulativeSlip, cumulative_slip, const SynchronizedArray &) AKANTU_GET_MACRO(SlipVelocity, slip_velocity, const SynchronizedArray &) /// set parameter of a given node /// (if you need to set to all: used the setMixed function of the Parsable). virtual void setParam(const std::string & name, UInt node, Real value); // replaced by the setMixed of the Parsable // virtual void setParam(const std::string & param, Real value) { // AKANTU_ERROR("Friction does not know the following parameter: " << // param); // }; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: - NTNBaseContact * contact; + NTNBaseContact & contact; // if node is sticking SynchronizedArray is_sticking; // frictional strength SynchronizedArray frictional_strength; // friction force SynchronizedArray friction_traction; // slip SynchronizedArray slip; SynchronizedArray cumulative_slip; // slip velocity (tangential vector) SynchronizedArray slip_velocity; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "ntn_base_friction_inline_impl.cc" /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const NTNBaseFriction & _this) { _this.printself(stream); return stream; } } // namespace akantu #endif /* __AST_NTN_BASE_FRICTION_HH__ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_contact.cc b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_contact.cc index 922b6f4c7..4e4232b23 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_contact.cc +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_contact.cc @@ -1,590 +1,554 @@ /** * @file ntn_contact.cc * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief implementation of ntn_contact * * @section LICENSE * * Copyright (©) 2015-2018 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 . * */ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_contact.hh" #include "dumper_nodal_field.hh" #include "dumper_text.hh" namespace akantu { /* -------------------------------------------------------------------------- */ -NTNContact::NTNContact(SolidMechanicsModel & model, const ContactID & id, +NTNContact::NTNContact(SolidMechanicsModel & model, const ID & id, const MemoryID & memory_id) : NTNBaseContact(model, id, memory_id), masters(0, 1, 0, id + ":masters", std::numeric_limits::quiet_NaN(), "masters"), lumped_boundary_masters(0, 1, 0, id + ":lumped_boundary_masters", std::numeric_limits::quiet_NaN(), "lumped_boundary_masters"), master_elements("master_elements", id, memory_id) { AKANTU_DEBUG_IN(); const Mesh & mesh = this->model.getMesh(); UInt spatial_dimension = this->model.getSpatialDimension(); - mesh.initElementTypeMapArray(this->master_elements, 1, spatial_dimension - 1); + this->master_elements.initialize(mesh, _nb_component = 1, + _spatial_dimension = spatial_dimension - 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::pairInterfaceNodes(const ElementGroup & slave_boundary, const ElementGroup & master_boundary, UInt surface_normal_dir, const Mesh & mesh, Array & pairs) { AKANTU_DEBUG_IN(); pairs.resize(0); AKANTU_DEBUG_ASSERT(pairs.getNbComponent() == 2, "Array of node pairs should have nb_component = 2," << " but has nb_component = " << pairs.getNbComponent()); UInt dim = mesh.getSpatialDimension(); AKANTU_DEBUG_ASSERT(surface_normal_dir < dim, "Mesh is of " << dim << " dimensions" << " and cannot have direction " << surface_normal_dir << " for surface normal"); // offset for projection computation - UInt offset[dim - 1]; + Vector offset(dim - 1); for (UInt i = 0, j = 0; i < dim; ++i) { if (surface_normal_dir != i) { - offset[j] = i; + offset(j) = i; ++j; } } // find projected node coordinates const Array & coordinates = mesh.getNodes(); // find slave nodes Array proj_slave_coord(slave_boundary.getNbNodes(), dim - 1, 0.); Array slave_nodes(slave_boundary.getNbNodes()); UInt n(0); - for (ElementGroup::const_node_iterator nodes_it(slave_boundary.node_begin()); - nodes_it != slave_boundary.node_end(); ++nodes_it) { + for (auto && slave_node : slave_boundary.getNodes()) { for (UInt d = 0; d < dim - 1; ++d) { - proj_slave_coord(n, d) = coordinates(*nodes_it, offset[d]); - slave_nodes(n) = *nodes_it; + proj_slave_coord(n, d) = coordinates(slave_node, offset[d]); + slave_nodes(n) = slave_node; } ++n; } // find master nodes Array proj_master_coord(master_boundary.getNbNodes(), dim - 1, 0.); Array master_nodes(master_boundary.getNbNodes()); n = 0; - for (ElementGroup::const_node_iterator nodes_it(master_boundary.node_begin()); - nodes_it != master_boundary.node_end(); ++nodes_it) { + for (auto && master_node : master_boundary.getNodes()) { for (UInt d = 0; d < dim - 1; ++d) { - proj_master_coord(n, d) = coordinates(*nodes_it, offset[d]); - master_nodes(n) = *nodes_it; + proj_master_coord(n, d) = coordinates(master_node, offset[d]); + master_nodes(n) = master_node; } ++n; } // find minimum distance between slave nodes to define tolerance Real min_dist = std::numeric_limits::max(); for (UInt i = 0; i < proj_slave_coord.size(); ++i) { for (UInt j = i + 1; j < proj_slave_coord.size(); ++j) { Real dist = 0.; for (UInt d = 0; d < dim - 1; ++d) { dist += (proj_slave_coord(i, d) - proj_slave_coord(j, d)) * (proj_slave_coord(i, d) - proj_slave_coord(j, d)); } if (dist < min_dist) { min_dist = dist; } } } min_dist = std::sqrt(min_dist); Real local_tol = 0.1 * min_dist; // find master slave node pairs for (UInt i = 0; i < proj_slave_coord.size(); ++i) { for (UInt j = 0; j < proj_master_coord.size(); ++j) { Real dist = 0.; for (UInt d = 0; d < dim - 1; ++d) { dist += (proj_slave_coord(i, d) - proj_master_coord(j, d)) * (proj_slave_coord(i, d) - proj_master_coord(j, d)); } dist = std::sqrt(dist); if (dist < local_tol) { // it is a pair - UInt pair[2]; + Vector pair(2); pair[0] = slave_nodes(i); pair[1] = master_nodes(j); pairs.push_back(pair); continue; // found master do not need to search further for this slave } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -void NTNContact::addSurfacePair(const Surface & slave, const Surface & master, +void NTNContact::addSurfacePair(const ID & slave, const ID & master, UInt surface_normal_dir) { AKANTU_DEBUG_IN(); const Mesh & mesh = this->model.getMesh(); const ElementGroup & slave_boundary = mesh.getElementGroup(slave); const ElementGroup & master_boundary = mesh.getElementGroup(master); this->contact_surfaces.insert(&slave_boundary); this->contact_surfaces.insert(&master_boundary); Array pairs(0, 2); NTNContact::pairInterfaceNodes(slave_boundary, master_boundary, surface_normal_dir, this->model.getMesh(), pairs); // eliminate pairs which contain a pbc slave node Array pairs_no_PBC_slaves(0, 2); Array::const_vector_iterator it = pairs.begin(2); Array::const_vector_iterator end = pairs.end(2); for (; it != end; ++it) { const Vector & pair = *it; - if (!this->model.isPBCSlaveNode(pair(0)) && - !this->model.isPBCSlaveNode(pair(1))) { + if (not mesh.isPeriodicSlave(pair(0)) and + not mesh.isPeriodicSlave(pair(1))) { pairs_no_PBC_slaves.push_back(pair); } } this->addNodePairs(pairs_no_PBC_slaves); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::addNodePairs(const Array & pairs) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(pairs.getNbComponent() == 2, "Array of node pairs should have nb_component = 2," << " but has nb_component = " << pairs.getNbComponent()); UInt nb_pairs = pairs.size(); for (UInt n = 0; n < nb_pairs; ++n) { this->addSplitNode(pairs(n, 0), pairs(n, 1)); } // synchronize with depending nodes findBoundaryElements(this->slaves.getArray(), this->slave_elements); findBoundaryElements(this->masters.getArray(), this->master_elements); updateInternalData(); syncArrays(_added); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::getNodePairs(Array & pairs) const { AKANTU_DEBUG_IN(); pairs.resize(0); AKANTU_DEBUG_ASSERT(pairs.getNbComponent() == 2, "Array of node pairs should have nb_component = 2," << " but has nb_component = " << pairs.getNbComponent()); UInt nb_pairs = this->getNbContactNodes(); for (UInt n = 0; n < nb_pairs; ++n) { - UInt pair[2]; - pair[0] = this->slaves(n); - pair[1] = this->masters(n); + Vector pair{this->slaves(n), this->masters(n)}; pairs.push_back(pair); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::addSplitNode(UInt slave, UInt master) { AKANTU_DEBUG_IN(); NTNBaseContact::addSplitNode(slave); this->masters.push_back(master); this->lumped_boundary_masters.push_back( std::numeric_limits::quiet_NaN()); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* This function only works for surface elements with one quad point. For surface elements with more quad points, it computes still, but the result might not be what you are looking for. */ void NTNContact::updateNormals() { AKANTU_DEBUG_IN(); // set normals to zero this->normals.clear(); // contact information UInt dim = this->model.getSpatialDimension(); UInt nb_contact_nodes = this->getNbContactNodes(); this->synch_registry->synchronize(_gst_cf_nodal); // synchronize current pos const Array & cur_pos = this->model.getCurrentPosition(); FEEngine & boundary_fem = this->model.getFEEngineBoundary(); const Mesh & mesh = this->model.getMesh(); - for (ghost_type_t::iterator gt = ghost_type_t::begin(); - gt != ghost_type_t::end(); ++gt) { - Mesh::type_iterator it = mesh.firstType(dim - 1, *gt); - Mesh::type_iterator last = mesh.lastType(dim - 1, *gt); - - for (; it != last; ++it) { + for (auto ghost_type: ghost_types) { + for (auto & type : mesh.elementTypes(dim - 1, ghost_type)) { // compute the normals Array quad_normals(0, dim); - boundary_fem.computeNormalsOnIntegrationPoints(cur_pos, quad_normals, *it, - *gt); + boundary_fem.computeNormalsOnIntegrationPoints(cur_pos, quad_normals, type, + ghost_type); - UInt nb_quad_points = boundary_fem.getNbIntegrationPoints(*it, *gt); + UInt nb_quad_points = boundary_fem.getNbIntegrationPoints(type, ghost_type); // new version: compute normals only based on master elements (and not all // boundary elements) // ------------------------------------------------------------------------------------- - UInt nb_nodes_per_element = mesh.getNbNodesPerElement(*it); - const Array & connectivity = mesh.getConnectivity(*it, *gt); + UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); + const Array & connectivity = mesh.getConnectivity(type, ghost_type); - Array::const_iterator elem_it = - (this->master_elements)(*it, *gt).begin(); - Array::const_iterator elem_it_end = - (this->master_elements)(*it, *gt).end(); // loop over contact nodes - for (; elem_it != elem_it_end; ++elem_it) { + for (auto & element : (this->master_elements)(type, ghost_type)) { for (UInt q = 0; q < nb_nodes_per_element; ++q) { - UInt node = connectivity(*elem_it, q); + UInt node = connectivity(element, q); UInt node_index = this->masters.find(node); - AKANTU_DEBUG_ASSERT(node_index != UInt(-1), - "Could not find node " << node - << " in the array!"); + AKANTU_DEBUG_ASSERT(node_index != UInt(-1), "Could not find node " + << node + << " in the array!"); for (UInt q = 0; q < nb_quad_points; ++q) { // add quad normal to master normal for (UInt d = 0; d < dim; ++d) { this->normals(node_index, d) += - quad_normals((*elem_it) * nb_quad_points + q, d); + quad_normals(element * nb_quad_points + q, d); } } } } - // ------------------------------------------------------------------------------------- - - /* - // get elements connected to each node - CSR node_to_element; - MeshUtils::buildNode2ElementsElementTypeMap(mesh, node_to_element, *it, - *gt); - - // add up normals to all master nodes - for (UInt n=0; nmasters(n); - CSR::iterator elem = node_to_element.begin(master); - // loop over all elements connected to this master node - for (; elem != node_to_element.end(master); ++elem) { - UInt e = *elem; - // loop over all quad points of this element - for (UInt q=0; qnormals(n,d) += quad_normals(e*nb_quad_points + q, d); - } - } - } - } - */ } } Real * master_normals = this->normals.storage(); for (UInt n = 0; n < nb_contact_nodes; ++n) { if (dim == 2) Math::normalize2(&(master_normals[n * dim])); else if (dim == 3) Math::normalize3(&(master_normals[n * dim])); } // // normalize normals // auto nit = this->normals.begin(); // auto nend = this->normals.end(); // for (; nit != nend; ++nit) { // nit->normalize(); // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::dumpRestart(const std::string & file_name) const { AKANTU_DEBUG_IN(); NTNBaseContact::dumpRestart(file_name); this->masters.dumpRestartFile(file_name); this->lumped_boundary_masters.dumpRestartFile(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::readRestart(const std::string & file_name) { AKANTU_DEBUG_IN(); NTNBaseContact::readRestart(file_name); this->masters.readRestartFile(file_name); this->lumped_boundary_masters.readRestartFile(file_name); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::updateImpedance() { AKANTU_DEBUG_IN(); UInt nb_contact_nodes = getNbContactNodes(); Real delta_t = this->model.getTimeStep(); AKANTU_DEBUG_ASSERT(delta_t != NAN, "Time step is NAN. Have you set it already?"); const Array & mass = this->model.getMass(); for (UInt n = 0; n < nb_contact_nodes; ++n) { UInt master = this->masters(n); UInt slave = this->slaves(n); Real imp = (this->lumped_boundary_masters(n) / mass(master)) + (this->lumped_boundary_slaves(n) / mass(slave)); imp = 2 / delta_t / imp; this->impedance(n) = imp; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::updateLumpedBoundary() { AKANTU_DEBUG_IN(); internalUpdateLumpedBoundary(this->slaves.getArray(), this->slave_elements, this->lumped_boundary_slaves); internalUpdateLumpedBoundary(this->masters.getArray(), this->master_elements, this->lumped_boundary_masters); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::applyContactPressure() { AKANTU_DEBUG_IN(); UInt nb_ntn_pairs = getNbContactNodes(); UInt dim = this->model.getSpatialDimension(); - Array & residual = this->model.getResidual(); + Array & residual = this->model.getInternalForce(); for (UInt n = 0; n < nb_ntn_pairs; ++n) { UInt master = this->masters(n); UInt slave = this->slaves(n); for (UInt d = 0; d < dim; ++d) { residual(master, d) += this->lumped_boundary_masters(n) * this->contact_pressure(n, d); residual(slave, d) -= this->lumped_boundary_slaves(n) * this->contact_pressure(n, d); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::computeRelativeTangentialField( const Array & field, Array & rel_tang_field) const { AKANTU_DEBUG_IN(); // resize arrays to zero rel_tang_field.resize(0); UInt dim = this->model.getSpatialDimension(); auto it_field = field.begin(dim); auto it_normal = this->normals.getArray().begin(dim); Vector rfv(dim); Vector np_rfv(dim); UInt nb_contact_nodes = this->slaves.size(); for (UInt n = 0; n < nb_contact_nodes; ++n) { // nodes UInt slave = this->slaves(n); UInt master = this->masters(n); // relative field vector (slave - master) rfv = Vector(it_field[slave]); rfv -= Vector(it_field[master]); // normal projection of relative field const Vector normal_v = it_normal[n]; np_rfv = normal_v; np_rfv *= rfv.dot(normal_v); // subract normal projection from relative field to get the tangential // projection rfv -= np_rfv; rel_tang_field.push_back(rfv); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::computeRelativeNormalField( const Array & field, Array & rel_normal_field) const { AKANTU_DEBUG_IN(); // resize arrays to zero rel_normal_field.resize(0); UInt dim = this->model.getSpatialDimension(); // Real * field_p = field.storage(); // Real * normals_p = this->normals.storage(); Array::const_iterator> it_field = field.begin(dim); Array::const_iterator> it_normal = this->normals.getArray().begin(dim); Vector rfv(dim); UInt nb_contact_nodes = this->getNbContactNodes(); for (UInt n = 0; n < nb_contact_nodes; ++n) { // nodes UInt slave = this->slaves(n); UInt master = this->masters(n); // relative field vector (slave - master) rfv = Vector(it_field[slave]); rfv -= Vector(it_field[master]); // length of normal projection of relative field const Vector normal_v = it_normal[n]; rel_normal_field.push_back(rfv.dot(normal_v)); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Int NTNContact::getNodeIndex(UInt node) const { AKANTU_DEBUG_IN(); Int slave_i = NTNBaseContact::getNodeIndex(node); Int master_i = this->masters.find(node); AKANTU_DEBUG_OUT(); return std::max(slave_i, master_i); } /* -------------------------------------------------------------------------- */ void NTNContact::printself(std::ostream & stream, int indent) const { AKANTU_DEBUG_IN(); std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "NTNContact [" << std::endl; NTNBaseContact::printself(stream, indent); stream << space << " + masters : " << std::endl; this->masters.printself(stream, indent + 2); stream << space << " + lumped_boundary_mastres : " << std::endl; this->lumped_boundary_masters.printself(stream, indent + 2); stream << space << "]" << std::endl; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::syncArrays(SyncChoice sync_choice) { AKANTU_DEBUG_IN(); NTNBaseContact::syncArrays(sync_choice); this->masters.syncElements(sync_choice); this->lumped_boundary_masters.syncElements(sync_choice); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NTNContact::addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id) { AKANTU_DEBUG_IN(); /* #ifdef AKANTU_USE_IOHELPER const Array & nodal_filter = this->slaves.getArray(); #define ADD_FIELD(field_id, field, type) \ internalAddDumpFieldToDumper(dumper_name, \ field_id, \ new DumperIOHelper::NodalField< type, true, \ Array, \ Array >(field, 0, 0, &nodal_filter)) */ if (field_id == "lumped_boundary_master") { internalAddDumpFieldToDumper( dumper_name, field_id, new dumper::NodalField(this->lumped_boundary_masters.getArray())); } else { NTNBaseContact::addDumpFieldToDumper(dumper_name, field_id); } /* #undef ADD_FIELD #endif */ AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_contact.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_contact.hh index ca9c16b90..2bc6d5a45 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_contact.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_contact.hh @@ -1,166 +1,166 @@ /** * @file ntn_contact.hh * * @author David Simon Kammer * * @date creation: Tue Dec 02 2014 * @date last modification: Fri Feb 23 2018 * * @brief contact for node to node discretization * * @section LICENSE * * Copyright (©) 2015-2018 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 __AST_NTN_CONTACT_HH__ #define __AST_NTN_CONTACT_HH__ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_base_contact.hh" namespace akantu { /* -------------------------------------------------------------------------- */ class NTNContact : public NTNBaseContact { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - NTNContact(SolidMechanicsModel & model, const ContactID & id = "contact", + NTNContact(SolidMechanicsModel & model, const ID & id = "contact", const MemoryID & memory_id = 0); virtual ~NTNContact(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// add surface pair and pair nodes according to the surface normal - void addSurfacePair(const Surface & slave, const Surface & master, + void addSurfacePair(const ID & slave, const ID & master, UInt surface_normal_dir); /// fills the pairs vector with interface node pairs (*,0)=slaves, /// (*,1)=masters static void pairInterfaceNodes(const ElementGroup & slave_boundary, const ElementGroup & master_boundary, UInt surface_normal_dir, const Mesh & mesh, Array & pairs); // add node pairs from a list with pairs(*,0)=slaves and pairs(*,1)=masters void addNodePairs(const Array & pairs); /// add node pair virtual void addSplitNode(UInt slave, UInt master); /// update (compute the normals on the master nodes) virtual void updateNormals(); /// update the lumped boundary B matrix virtual void updateLumpedBoundary(); /// update the impedance matrix virtual void updateImpedance(); /// impose the normal contact force virtual void applyContactPressure(); /// dump restart file virtual void dumpRestart(const std::string & file_name) const; /// read restart file virtual void readRestart(const std::string & file_name); /// compute the normal gap virtual void computeNormalGap(Array & gap) const { this->computeRelativeNormalField(this->model.getCurrentPosition(), gap); }; /// compute relative normal field (only value that has to be multiplied with /// the normal) /// relative to master nodes virtual void computeRelativeNormalField(const Array & field, Array & rel_normal_field) const; /// compute relative tangential field (complet array) /// relative to master nodes virtual void computeRelativeTangentialField(const Array & field, Array & rel_tang_field) const; /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; protected: /// synchronize arrays virtual void syncArrays(SyncChoice sync_choice); /* ------------------------------------------------------------------------ */ /* Dumpable */ /* ------------------------------------------------------------------------ */ public: virtual void addDumpFieldToDumper(const std::string & dumper_name, const std::string & field_id); // virtual void addDumpFieldVector(const std::string & field_id); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Masters, masters, const SynchronizedArray &) AKANTU_GET_MACRO(LumpedBoundaryMasters, lumped_boundary_masters, const SynchronizedArray &) /// get interface node pairs (*,0) are slaves, (*,1) are masters void getNodePairs(Array & pairs) const; /// get index of node in either slaves or masters array /// if node is in neither of them, return -1 Int getNodeIndex(UInt node) const; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// array of master nodes SynchronizedArray masters; /// lumped boundary of master nodes SynchronizedArray lumped_boundary_masters; // element list for dump and lumped_boundary ElementTypeMapArray master_elements; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "ntn_contact_inline_impl.cc" /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const NTNContact & _this) { _this.printself(stream); return stream; } } // namespace akantu #endif /* __AST_NTN_CONTACT_HH__ */ diff --git a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_friction.hh b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_friction.hh index 669c13c2c..149295355 100644 --- a/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_friction.hh +++ b/extra_packages/traction-at-split-node-contact/src/ntn_contact/ntn_friction.hh @@ -1,100 +1,100 @@ /** * @file ntn_friction.hh * * @author David Simon Kammer * * @date creation: Fri Jun 18 2010 * @date last modification: Fri Feb 23 2018 * * @brief implementation of friction for node to node contact * * @section LICENSE * * Copyright (©) 2015-2018 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 __AST_NTN_FRICTION_HH__ #define __AST_NTN_FRICTION_HH__ /* -------------------------------------------------------------------------- */ // simtools #include "ntn_base_friction.hh" #include "ntn_friclaw_coulomb.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template