diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt
index 04b727722..c7274da8d 100644
--- a/python/CMakeLists.txt
+++ b/python/CMakeLists.txt
@@ -1,301 +1,305 @@
#===============================================================================
# @file CMakeLists.txt
#
# @author Nicolas Richart <nicolas.richart@epfl.ch>
#
# @date creation: Fri Dec 12 2014
# @date last modification: Mon Jan 18 2016
#
# @brief CMake file for the python wrapping of akantu
#
# @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 <http://www.gnu.org/licenses/>.
#
#===============================================================================
#===============================================================================
# Configuration
#===============================================================================
package_get_all_definitions(AKANTU_DEFS)
list(REMOVE_ITEM AKANTU_DEFS AKANTU_CORE_CXX11)
#message(${AKANTU_DEFS})
set(AKA_DEFS "")
foreach (def ${AKANTU_DEFS})
list(APPEND AKA_DEFS "-D${def}")
endforeach()
set(AKANTU_SWIG_FLAGS -w309,325,401,317,509,503,383,384 ${AKA_DEFS})
set(AKANTU_SWIG_OUTDIR ${CMAKE_CURRENT_SOURCE_DIR})
set(AKANTU_SWIG_MODULES swig/akantu.i)
#===============================================================================
# Swig wrapper
#===============================================================================
set(SWIG_REQURIED_VERISON 3.0)
find_package(SWIG ${SWIG_REQURIED_VERISON})
find_package(PythonInterp 2.7 REQUIRED)
mark_as_advanced(SWIG_EXECUTABLE)
package_get_all_include_directories(
AKANTU_LIBRARY_INCLUDE_DIRS
)
package_get_all_external_informations(
AKANTU_EXTERNAL_INCLUDE_DIR
AKANTU_EXTERNAL_LIBRARIES
)
set(_swig_include_dirs
${CMAKE_CURRENT_SOURCE_DIR}/swig
${AKANTU_LIBRARY_INCLUDE_DIRS}
${PROJECT_BINARY_DIR}/src
${AKANTU_EXTERNAL_INCLUDE_DIR}
)
include(CMakeParseArguments)
function(swig_generate_dependencies _module _depedencies _depedencies_out)
set(_dependencies_script "${CMAKE_CURRENT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/_swig_generate_dependencies.cmake")
file(WRITE ${_dependencies_script} "
set(_include_directories ${_include_directories})
list(APPEND _include_directories \"./\")
set(_dep)
set(_files_to_process \${_module})
while(_files_to_process)
list(GET _files_to_process 0 _file)
list(REMOVE_AT _files_to_process 0)
file(STRINGS \${_file} _file_content REGEX \"^%include *\\\"(.*)\\\"\")
set(_includes)
foreach(_line \${_file_content})
string(REGEX REPLACE \"^%include *\\\"(.*)\\\"\" \"\\\\1\" _inc \${_line})
if(_inc)
list(APPEND _includes \${_inc})
endif()
endforeach()
foreach(_include \${_includes})
unset(_found)
foreach(_inc_dir \${_include_directories})
if(EXISTS \${_inc_dir}/\${_include})
set(_found \${_inc_dir}/\${_include})
break()
endif()
endforeach()
if(_found)
list(APPEND _files_to_process \${_found})
list(APPEND _dep \${_found})
endif()
endforeach()
endwhile()
get_filename_component(_module_we \"\${_module}\" NAME_WE)
set(_dependencies_file \${CMAKE_CURRENT_BINARY_DIR}\${CMAKE_FILES_DIRECTORY}/_swig_\${_module_we}_depends.cmake)
file(WRITE \"\${_dependencies_file}\"
\"set(_swig_\${_module_we}_depends\")
foreach(_d \${_dep})
file(APPEND \"\${_dependencies_file}\" \"
\${_d}\")
endforeach()
file(APPEND \"\${_dependencies_file}\" \"
)\")
")
get_filename_component(_module_absolute "${_module}" ABSOLUTE)
get_filename_component(_module_we "${_module}" NAME_WE)
set(_dependencies_file ${CMAKE_CURRENT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/_swig_${_module_we}_depends.cmake)
if(EXISTS ${_dependencies_file})
include(${_dependencies_file})
else()
execute_process(COMMAND ${CMAKE_COMMAND}
-D_module=${_module_absolute}
-P ${_dependencies_script}
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
include(${_dependencies_file})
endif()
set(${_depedencies_out} ${_swig_${_module_we}_depends} PARENT_SCOPE)
add_custom_command(OUTPUT ${_dependencies_file}
COMMAND ${CMAKE_COMMAND}
-D_module=${_module_absolute}
-P ${_dependencies_script}
COMMENT "Scanning dependencies for swig module ${_module_we}"
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
MAIN_DEPENDENCY ${_module_absolute}
DEPENDS ${_swig_${_module_we}_depends}
)
set(${_depedencies} ${_dependencies_file} PARENT_SCOPE)
endfunction()
function(swig_generate_wrappers project _wrappers_cpp _wrappers_py)
cmake_parse_arguments(_swig_opt "" "OUTPUT_DIR;DEPENDENCIES" "EXTRA_FLAGS;INCLUDE_DIRECTORIES" ${ARGN})
if(_swig_opt_OUTPUT_DIR)
set(_output_dir ${_swig_opt_OUTPUT_DIR})
else()
set(_output_dir ${CMAKE_CURRENT_BINARY_DIR})
endif()
set(_swig_wrappers)
get_directory_property(_include_directories INCLUDE_DIRECTORIES)
list(APPEND _include_directories ${_swig_opt_INCLUDE_DIRECTORIES})
if(_include_directories)
string(REPLACE ";" ";-I" _swig_include_directories "${_include_directories}")
endif()
foreach(_module ${_swig_opt_UNPARSED_ARGUMENTS})
swig_generate_dependencies(${_module} _module_dependencies _depends_out)
if(_swig_opt_DEPENDENCIES)
set(${_swig_opt_DEPENDENCIES} ${_depends_out} PARENT_SCOPE)
endif()
get_filename_component(_module_absolute "${_module}" ABSOLUTE)
get_filename_component(_module_path "${_module_absolute}" PATH)
get_filename_component(_module_name "${_module}" NAME)
get_filename_component(_module_we "${_module}" NAME_WE)
set(_wrapper "${_output_dir}/${_module_we}_wrapper.cpp")
set(_extra_wrapper "${_output_dir}/${_module_we}.py")
set(_extra_wrapper_bin "${CMAKE_CURRENT_BINARY_DIR}/${_module_we}.py")
if(SWIG_FOUND)
set_source_files_properties("${_wrapper}" PROPERTIES GENERATED 1)
set_source_files_properties("${_extra_wrapper}" PROPERTIES GENERATED 1)
set(_dependencies_file ${CMAKE_CURRENT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/_swig_${_module_we}_depends.cmake)
set(_ouput "${_wrapper}" "${_extra_wrapper}")
add_custom_command(
OUTPUT ${_ouput}
COMMAND "${SWIG_EXECUTABLE}"
ARGS -python -c++
${_swig_opt_EXTRA_FLAGS}
-outdir ${_output_dir}
-I${_swig_include_directories} -I${_module_path}
-o "${_wrapper}"
"${_module_absolute}"
COMMAND ${CMAKE_COMMAND} -E copy_if_different ${_extra_wrapper} ${_extra_wrapper_bin}
# MAIN_DEPENDENCY "${_module_absolute}"
DEPENDS ${_module_dependencies}
COMMENT "Generating swig wrapper ${_module} -> ${_wrapper}"
)
list(APPEND _swig_wrappers ${_wrapper})
list(APPEND _swig_wrappers_py "${_extra_wrapper_bin}")
else()
if(NOT EXISTS ${_wrapper} OR NOT EXISTS "${_extra_wrapper}")
message(FATAL_ERROR "The file ${_wrapper} and/or ${_extra_wrapper} does "
"not exists and they cannot be generated. Install swig ${SWIG_REQURIED_VERISON} "
" in order to generate them. Or get them from a different machine, "
"in order to be able to compile the python interface")
else()
list(APPEND _swig_wrappers "${_wrapper}")
list(APPEND _swig_wrappers_py "${_extra_wrapper_bin}")
endif()
endif()
endforeach()
add_custom_target(${project}_generate_swig_wrappers DEPENDS ${_swig_wrappers})
set(${_wrappers_cpp} ${_swig_wrappers} PARENT_SCOPE)
set(${_wrappers_py} ${_swig_wrappers_py} PARENT_SCOPE)
endfunction()
swig_generate_wrappers(akantu AKANTU_SWIG_WRAPPERS_CPP AKANTU_WRAPPERS_PYTHON
${AKANTU_SWIG_MODULES}
EXTRA_FLAGS ${AKANTU_SWIG_FLAGS}
DEPENDENCIES _deps
INCLUDE_DIRECTORIES ${_swig_include_dirs})
if(AKANTU_SWIG_WRAPPERS_CPP)
string(REPLACE ";" "', '" _ext_files "${AKANTU_SWIG_WRAPPERS_CPP}")
string(REPLACE ";" "', '" _inc_dirs "${_swig_include_dirs}")
- string(STRIP "${AKANTU_EXTRA_CXX_FLAGS}" _tmp_flags)
- string(REGEX REPLACE " +" "', '" _flags "${_tmp_flags} -Wno-maybe-uninitialized")
+ if (AKANTU_EXTRA_CXX_FLAGS)
+ string(STRIP "${AKANTU_EXTRA_CXX_FLAGS}" _tmp_flags)
+ string(REGEX REPLACE " +" "', '" _flags "${_tmp_flags} -Wno-maybe-uninitialized")
+ else()
+ set(_flags "-Wno-maybe-uninitialized")
+ endif()
if(CMAKE_VERBOSE_MAKEFILE)
set(_quiet)
else()
set(_quiet --quiet)
endif()
file(WRITE ${CMAKE_CURRENT_BINARY_DIR}/setup.py "
from distutils.core import setup
from distutils.core import setup, Extension
import os
os.environ['CC'] = '${CMAKE_CXX_COMPILER}'
os.environ['CXX'] = '${CMAKE_CXX_COMPILER}'
setup(name='akantu',
license='LGPLv3',
version='${AKANTU_VERSION}',
py_modules=['akantu'],
ext_modules=[Extension('_akantu', ['${_ext_files}'],
include_dirs=['${_inc_dirs}'],
language='c++',
libraries=['akantu'],
library_dirs=['${PROJECT_BINARY_DIR}/src'],
- extra_compile_args=['${_flags}'])],
+ extra_compile_args=['${_flags}', '-std=c++14'])],
)
")
add_custom_target(_akantu ALL
COMMAND ${PYTHON_EXECUTABLE} ./setup.py ${_quiet} --no-user-cfg build_ext --inplace
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
DEPENDS ${AKANTU_SWIG_WRAPPERS_CPP} akantu
COMMENT "Building akantu's python interface"
)
set_directory_properties(PROPERTIES
ADDITIONAL_MAKE_CLEAN_FILES
${CMAKE_CURRENT_BINARY_DIR}/_akantu${CMAKE_SHARED_MODULE_SUFFIX})
# add_library(_akantu MODULE ${AKANTU_SWIG_WRAPPERS_CPP})
# target_link_libraries(_akantu akantu)
# set_target_properties(_akantu PROPERTIES PREFIX "")
# list(APPEND AKANTU_EXPORT_LIST _akantu)
# install(TARGETS _akantu
# EXPORT ${AKANTU_TARGETS_EXPORT}
# LIBRARY DESTINATION lib COMPONENT python NAMELINK_SKIP # for real systems
# ARCHIVE DESTINATION lib COMPONENT python
# RUNTIME DESTINATION bin COMPONENT python # for windows ...
# )
# install(FILES ${AKANTU_WRAPPERS_PYTHON}
# DESTINATION lib COMPONENT python
# )
install(CODE "execute_process(
COMMAND ${PYTHON_EXECUTABLE} ./setup.py ${_quiet} install --prefix=${AKANTU_PYTHON_INSTALL_PREFIX}
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
)"
COMPONENT python)
package_declare_extra_files_to_package(python_interface ${_deps} ${PROJECT_SOURCE_DIR}/python/${AKANTU_SWIG_MODULES})
endif()
diff --git a/python/swig/aka_array.i b/python/swig/aka_array.i
index 117f261fb..5012a68a5 100644
--- a/python/swig/aka_array.i
+++ b/python/swig/aka_array.i
@@ -1,248 +1,251 @@
/**
* @file aka_array.i
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Fri Dec 12 2014
* @date last modification: Wed Nov 11 2015
*
* @brief wrapper for arrays
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
%{
#define SWIG_FILE_WITH_INIT
#include "aka_array.hh"
#include "aka_types.hh"
%}
%include "typemaps.i"
namespace akantu {
%ignore Array::operator=;
%ignore Array::operator[];
%ignore Array::operator();
%ignore Array::set;
%ignore Array::begin;
%ignore Array::end;
%ignore Array::begin_reinterpret;
%ignore Array::end_reinterpret;
+ %ignore ArrayBase::getSize;
+ %ignore Array::operator*=;
+ %ignore Array::operator*;
};
%include "aka_array.hh"
namespace akantu {
%ignore TensorProxy::operator=;
%ignore TensorProxy::operator[];
%ignore TensorProxy::operator();
%ignore Tensor3Proxy::operator=;
%ignore Tensor3Proxy::operator[];
%ignore Tensor3Proxy::operator();
%ignore TensorStorage::operator=;
%ignore TensorStorage::operator[];
%ignore TensorStorage::operator();
%ignore VectorProxy::operator=;
%ignore VectorProxy::operator[];
%ignore VectorProxy::operator();
%ignore MatrixProxy::operator=;
%ignore MatrixProxy::operator[];
%ignore MatrixProxy::operator();
%ignore Matrix::operator=;
%ignore Matrix::operator[];
%ignore Matrix::operator();
%ignore Tensor3::operator=;
%ignore Tensor3::operator[];
%ignore Tensor3::operator();
%ignore Vector::operator=;
%ignore Vector::operator[];
%ignore Vector::operator();
%ignore Vector::solve;
};
%include "aka_types.hh"
namespace akantu {
%template(RArray) Array<akantu::Real, true>;
%template(UArray) Array<akantu::UInt, true>;
%template(BArray) Array<bool, true>;
%template(RVector) Vector<akantu::Real>;
};
%include "numpy.i"
%init %{
import_array();
%}
%inline %{
namespace akantu{
template <typename T>
class ArrayForPython : public Array<T>{
public:
ArrayForPython(T * wrapped_memory,
- UInt size = 0,
+ UInt size_ = 0,
UInt nb_component = 1,
const ID & id = "")
: Array<T>(0,nb_component,id){
this->values = wrapped_memory;
- this->size = size;
+ this->size_ = size_;
};
~ArrayForPython(){
this->values = NULL;
};
void resize(UInt new_size){
- AKANTU_DEBUG_ASSERT(this->size == new_size,"cannot resize a temporary vector");
+ AKANTU_DEBUG_ASSERT(this->size_ == new_size,"cannot resize a temporary vector");
}
};
}
template <typename T> int getPythonDataTypeCode(){ AKANTU_EXCEPTION("undefined type");}
template <> int getPythonDataTypeCode<bool>(){
int data_typecode = NPY_NOTYPE;
size_t s = sizeof(bool);
switch(s) {
case 1: data_typecode = NPY_BOOL; break;
case 2: data_typecode = NPY_UINT16; break;
case 4: data_typecode = NPY_UINT32; break;
case 8: data_typecode = NPY_UINT64; break;
}
return data_typecode;
}
template <> int getPythonDataTypeCode<double>(){return NPY_DOUBLE;}
template <> int getPythonDataTypeCode<long double>(){return NPY_LONGDOUBLE;}
template <> int getPythonDataTypeCode<float>(){return NPY_FLOAT;}
template <> int getPythonDataTypeCode<unsigned long>(){
int data_typecode = NPY_NOTYPE;
size_t s = sizeof(unsigned long);
switch(s) {
case 2: data_typecode = NPY_UINT16; break;
case 4: data_typecode = NPY_UINT32; break;
case 8: data_typecode = NPY_UINT64; break;
}
return data_typecode;
}
template <> int getPythonDataTypeCode<akantu::UInt>(){
int data_typecode = NPY_NOTYPE;
size_t s = sizeof(akantu::UInt);
switch(s) {
case 2: data_typecode = NPY_UINT16; break;
case 4: data_typecode = NPY_UINT32; break;
case 8: data_typecode = NPY_UINT64; break;
}
return data_typecode;
}
template <> int getPythonDataTypeCode<int>(){
int data_typecode = NPY_NOTYPE;
size_t s = sizeof(int);
switch(s) {
case 2: data_typecode = NPY_INT16; break;
case 4: data_typecode = NPY_INT32; break;
case 8: data_typecode = NPY_INT64; break;
}
return data_typecode;
}
int getSizeOfPythonType(int type_num){
switch (type_num){
case NPY_INT16 : return 2;break;
case NPY_UINT16: return 2;break;
case NPY_INT32 : return 4;break;
case NPY_UINT32: return 4;break;
case NPY_INT64 : return 8;break;
case NPY_UINT64: return 8;break;
case NPY_FLOAT: return sizeof(float);break;
case NPY_DOUBLE: return sizeof(double);break;
case NPY_LONGDOUBLE: return sizeof(long double);break;
}
return 0;
}
std::string getPythonTypeName(int type_num){
switch (type_num){
case NPY_INT16 : return "NPY_INT16" ;break;
case NPY_UINT16: return "NPY_UINT16";break;
case NPY_INT32 : return "NPY_INT32" ;break;
case NPY_UINT32: return "NPY_UINT32";break;
case NPY_INT64 : return "NPY_INT64" ;break;
case NPY_UINT64: return "NPY_UINT64";break;
case NPY_FLOAT: return "NPY_FLOAT" ;break;
case NPY_DOUBLE: return "NPY_DOUBLE";break;
case NPY_LONGDOUBLE: return "NPY_LONGDOUBLE";break;
}
return 0;
}
template <typename T> void checkDataType(int type_num){
AKANTU_DEBUG_ASSERT(type_num == getPythonDataTypeCode<T>(),
"incompatible types between numpy and input function: " <<
type_num << " != " << getPythonDataTypeCode<T>() << std::endl <<
getSizeOfPythonType(type_num) << " != " << sizeof(T) <<
std::endl <<
"The numpy array is of type " << getPythonTypeName(type_num));
}
%}
%define %akantu_array_typemaps(DATA_TYPE)
%typemap(out, fragment="NumPy_Fragments") (akantu::Array< DATA_TYPE > &)
{
int data_typecode = getPythonDataTypeCode< DATA_TYPE >();
- npy_intp dims[2] = {npy_intp($1->getSize()), npy_intp($1->getNbComponent())};
+ npy_intp dims[2] = {npy_intp($1->size()), npy_intp($1->getNbComponent())};
PyObject* obj = PyArray_SimpleNewFromData(2, dims, data_typecode, $1->storage());
PyArrayObject* array = (PyArrayObject*) obj;
if (!array) SWIG_fail;
$result = SWIG_Python_AppendOutput($result, obj);
}
%typemap(in) akantu::Array< DATA_TYPE > &
{
if (!PyArray_Check($input)) {
AKANTU_EXCEPTION("incompatible input which is not a numpy");
}
else {
PyArray_Descr * numpy_type = (PyArray_Descr*)PyArray_DESCR((PyArrayObject*)$input);
int type_num = numpy_type->type_num;
checkDataType< DATA_TYPE >(type_num);
UInt _n = PyArray_NDIM((PyArrayObject*)$input);
if (_n != 2) AKANTU_EXCEPTION("incompatible numpy dimension " << _n);
npy_intp * ndims = PyArray_DIMS((PyArrayObject*)$input);
akantu::UInt sz = ndims[0];
akantu::UInt nb_components = ndims[1];
PyArrayIterObject *iter = (PyArrayIterObject *)PyArray_IterNew($input);
if (iter == NULL) AKANTU_EXCEPTION("Python internal error");
$1 = new akantu::ArrayForPython< DATA_TYPE >((DATA_TYPE*)(iter->dataptr),sz,nb_components,"tmp_array_for_python");
}
}
%enddef
%akantu_array_typemaps(double )
%akantu_array_typemaps(float )
%akantu_array_typemaps(unsigned int)
%akantu_array_typemaps(unsigned long)
%akantu_array_typemaps(int )
%akantu_array_typemaps(bool )
diff --git a/python/swig/aka_common.i b/python/swig/aka_common.i
index 95984de36..239a2fbc6 100644
--- a/python/swig/aka_common.i
+++ b/python/swig/aka_common.i
@@ -1,124 +1,125 @@
/**
* @file aka_common.i
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Fabian Barras <fabian.barras@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Fri Dec 12 2014
* @date last modification: Wed Jan 13 2016
*
* @brief wrapper to aka_common.hh
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
%{
- #include "aka_common.hh"
- #include "aka_csr.hh"
- #include "element.hh"
+#include "aka_common.hh"
+#include "aka_csr.hh"
+#include "element.hh"
%}
namespace akantu {
%ignore getStaticParser;
%ignore getUserParser;
+ %ignore ghost_types;
%ignore initialize(int & argc, char ** & argv);
%ignore initialize(const std::string & input_file, int & argc, char ** & argv);
extern const Array<UInt> empty_filter;
}
%typemap(in) (int argc, char *argv[]) {
int i = 0;
if (!PyList_Check($input)) {
PyErr_SetString(PyExc_ValueError, "Expecting a list");
return NULL;
}
$1 = PyList_Size($input);
$2 = new char *[$1+1];
for (i = 0; i < $1; i++) {
PyObject *s = PyList_GetItem($input,i);
if (!PyString_Check(s)) {
free($2);
PyErr_SetString(PyExc_ValueError, "List items must be strings");
return NULL;
}
$2[i] = PyString_AsString(s);
}
$2[i] = 0;
}
%typemap(freearg) (int argc, char *argv[]) {
%#if defined(__INTEL_COMPILER)
//#pragma warning ( disable : 383 )
%#elif defined (__clang__) // test clang to be sure that when we test for gnu it is only gnu
%#elif (defined(__GNUC__) || defined(__GNUG__))
%# if __cplusplus > 199711L
%# pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
%# endif
%#endif
delete [] $2;
%#if defined(__INTEL_COMPILER)
//#pragma warning ( disable : 383 )
%#elif defined (__clang__) // test clang to be sure that when we test for gnu it is only gnu
%#elif (defined(__GNUC__) || defined(__GNUG__))
%# if __cplusplus > 199711L
%# pragma GCC diagnostic pop
%# endif
%#endif
}
%inline %{
namespace akantu {
#if defined(AKANTU_USE_MPI)
const int MPI=1;
#else
const int MPI=0;
#endif
void _initializeWithArgv(const std::string & input_file, int argc, char *argv[]) {
initialize(input_file, argc, argv);
}
}
%}
%pythoncode %{
import sys as _aka_sys
def initialize(input_file="", argv=_aka_sys.argv):
if _aka_sys.modules[__name__].MPI == 1:
try:
from mpi4py import MPI
except ImportError:
pass
_initializeWithArgv(input_file, argv)
%}
%include "aka_config.hh"
%include "aka_common.hh"
%include "aka_element_classes_info.hh"
%include "element.hh"
diff --git a/python/swig/material.i b/python/swig/material.i
index 904b44bfa..c1d3b152c 100644
--- a/python/swig/material.i
+++ b/python/swig/material.i
@@ -1,49 +1,50 @@
/**
* @file material.i
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @brief material wrapper
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
%{
#include "solid_mechanics_model.hh"
#include "material_python.hh"
%}
namespace akantu {
%ignore Material::onNodesAdded;
%ignore Material::onNodesRemoved;
%ignore Material::onElementsAdded;
%ignore Material::onElementsRemoved;
%ignore Material::onElementsChanged;
+ %ignore Material::getParam;
}
%include "material.hh"
%extend akantu::Material {
Array<Real> & getArrayReal(const ID & id, const ElementType & type,
const GhostType & ghost_type = _not_ghost) {
return $self->getArray<Real>(id, type, ghost_type);
}
}
diff --git a/python/swig/mesh.i b/python/swig/mesh.i
index 57ede82de..53b7e2561 100644
--- a/python/swig/mesh.i
+++ b/python/swig/mesh.i
@@ -1,177 +1,177 @@
/**
* @file mesh.i
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Fabian Barras <fabian.barras@epfl.ch>
* @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Fri Dec 12 2014
* @date last modification: Wed Jan 13 2016
*
* @brief mesh wrapper
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
%{
#include "mesh.hh"
#include "node_group.hh"
#include "solid_mechanics_model.hh"
#include "dumpable_inline_impl.hh"
using akantu::IntegrationPoint;
using akantu::Vector;
using akantu::ElementTypeMapArray;
using akantu::MatrixProxy;
using akantu::Matrix;
using akantu::UInt;
using akantu::Real;
using akantu::Array;
using akantu::SolidMechanicsModel;
%}
namespace akantu {
%ignore NewNodesEvent;
%ignore RemovedNodesEvent;
%ignore NewElementsEvent;
%ignore RemovedElementsEvent;
%ignore MeshEventHandler;
%ignore MeshEvent< UInt >;
%ignore MeshEvent< Element >;
%ignore Mesh::extractNodalCoordinatesFromPBCElement;
%ignore Mesh::getGroupDumer;
%ignore Mesh::getFacetLocalConnectivity;
%ignore Mesh::getAllFacetTypes;
%ignore Mesh::getCommunicator;
}
print_self(Mesh)
// Swig considers enums to be ints, and it creates a conflict with two versions of getNbElement()
%rename(getNbElementByDimension)
akantu::Mesh::getNbElement(const UInt spatial_dimension = _all_dimensions,
const GhostType& ghost_type = _not_ghost,
const ElementKind& kind = _ek_not_defined) const;
%extend akantu::Mesh {
void resizeMesh(UInt nb_nodes, UInt nb_element, const ElementType & type) {
Array<Real> & nodes = const_cast<Array<Real> &>($self->getNodes());
nodes.resize(nb_nodes);
$self->addConnectivityType(type);
Array<UInt> & connectivity = const_cast<Array<UInt> &>($self->getConnectivity(type));
connectivity.resize(nb_element);
}
#if defined(AKANTU_COHESIVE_ELEMENT)
Array<Real> & getCohesiveBarycenter(SpacialDirection dir) {
UInt spatial_dimension = $self->getSpatialDimension();
ElementTypeMapArray<Real> & barycenter =
$self->registerData<Real>("barycenter");
$self->initElementTypeMapArray(barycenter, 1, spatial_dimension, false,
akantu::_ek_cohesive, true);
akantu::ElementType type = *($self->firstType(
spatial_dimension, akantu::_not_ghost, akantu::_ek_cohesive));
Vector<Real> bary(spatial_dimension);
Array<Real> & bary_coh = barycenter(type);
for (UInt i = 0; i < $self->getNbElement(type); ++i) {
bary.clear();
$self->getBarycenter(i, type, bary.storage());
bary_coh(i) = bary(dir);
}
return bary_coh;
}
#endif
}
%extend akantu::GroupManager {
void createGroupsFromStringMeshData(const std::string & dataset_name) {
$self->createGroupsFromMeshData<std::string>(dataset_name);
}
void createGroupsFromUIntMeshData(const std::string & dataset_name) {
$self->createGroupsFromMeshData<akantu::UInt>(dataset_name);
}
}
%extend akantu::NodeGroup {
akantu::Array<akantu::Real> & getGroupedNodes(akantu::Array<akantu::Real, true> & surface_array, Mesh & mesh) {
akantu::Array<akantu::UInt> group_node = $self->getNodes();
akantu::Array<akantu::Real> & full_array = mesh.getNodes();
- surface_array.resize(group_node.getSize());
+ surface_array.resize(group_node.size());
- for (UInt i = 0; i < group_node.getSize(); ++i) {
+ for (UInt i = 0; i < group_node.size(); ++i) {
for (UInt cmp = 0; cmp < full_array.getNbComponent(); ++cmp) {
surface_array(i,cmp) = full_array(group_node(i),cmp);
}
}
akantu::Array<akantu::Real> & res(surface_array);
return res;
}
akantu::Array<akantu::Real> & getGroupedArray(akantu::Array<akantu::Real, true> & surface_array, akantu::SolidMechanicsModel & model, int type) {
akantu::Array<akantu::Real> * full_array;
switch (type) {
case 0 : full_array = new akantu::Array<akantu::Real>(model.getDisplacement());
break;
case 1 : full_array = new akantu::Array<akantu::Real>(model.getVelocity());
break;
case 2 : full_array = new akantu::Array<akantu::Real>(model.getForce());
break;
}
akantu::Array<akantu::UInt> group_node = $self->getNodes();
- surface_array.resize(group_node.getSize());
+ surface_array.resize(group_node.size());
- for (UInt i = 0; i < group_node.getSize(); ++i) {
+ for (UInt i = 0; i < group_node.size(); ++i) {
for (UInt cmp = 0; cmp < full_array->getNbComponent(); ++cmp) {
surface_array(i,cmp) = (*full_array)(group_node(i),cmp);
}
}
akantu::Array<akantu::Real> & res(surface_array);
return res;
}
}
%include "group_manager.hh"
%include "element_group.hh"
%include "node_group.hh"
%include "dumper_iohelper.hh"
%include "dumpable_iohelper.hh"
%include "mesh.hh"
namespace akantu{
%extend Dumpable {
void addDumpFieldExternalReal(const std::string & field_id,
const Array<Real> & field){
$self->addDumpFieldExternal<Real>(field_id,field);
}
}
}
diff --git a/python/swig/mesh_utils.i b/python/swig/mesh_utils.i
index e5984ff34..31a551364 100644
--- a/python/swig/mesh_utils.i
+++ b/python/swig/mesh_utils.i
@@ -1,40 +1,43 @@
/**
* @file mesh_utils.i
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Fri Dec 12 2014
* @date last modification: Thu Jul 23 2015
*
* @brief mesh_utils wrapper
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
-namespace akantu {
- %ignore MeshPartition::getPartitions;
- %ignore MeshPartition::getPartition;
- %ignore MeshPartition::getGhostPartitionCSR;
-}
+%{
+#include "mesh_utils.hh"
+%}
+/* namespace akantu { */
+/* %ignore MeshPartition::getPartitions; */
+/* %ignore MeshPartition::getPartition; */
+/* %ignore MeshPartition::getGhostPartitionCSR; */
+/* } */
-%include "mesh_partition.hh"
+/* %include "mesh_partition.hh" */
%include "mesh_utils.hh"
diff --git a/python/swig/model.i b/python/swig/model.i
index 8775f3504..d08eaa605 100644
--- a/python/swig/model.i
+++ b/python/swig/model.i
@@ -1,80 +1,85 @@
/**
* @file model.i
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Fri Dec 12 2014
* @date last modification: Wed Nov 11 2015
*
* @brief model wrapper
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
%{
#include "boundary_condition_python_functor.hh"
%}
namespace akantu {
%ignore Model::createSynchronizerRegistry;
+ %ignore Model::getSynchronizerRegistry;
%ignore Model::createParallelSynch;
%ignore Model::getDOFSynchronizer;
- //%ignore Model::getSynchronizerRegistry;
%ignore Model::registerFEEngineObject;
%ignore Model::unregisterFEEngineObject;
%ignore Model::getFEEngineBoundary;
// %ignore Model::getFEEngine;
%ignore Model::getFEEngineClass;
%ignore Model::getFEEngineClassBoundary;
%ignore Model::setParser;
-
+ %ignore Model::updateDataForNonLocalCriterion;
%ignore IntegrationPoint::operator=;
%ignore FEEngine::getNbIntegrationPoints;
%ignore FEEngine::getShapes;
%ignore FEEngine::getShapesDerivatives;
%ignore FEEngine::getIntegrationPoints;
%ignore FEEngine::getIGFEMElementTypes;
%ignore FEEngine::interpolateOnIntegrationPoints(const Array<Real> &,ElementTypeMapArray<Real> &,const ElementTypeMapArray<UInt> *) const;
%ignore FEEngine::interpolateOnIntegrationPoints(const Array<Real> &,ElementTypeMapArray<Real> &) const;
%ignore FEEngine::interpolateOnIntegrationPoints(const Array<Real> &,Array<Real> &,UInt,const ElementType&,const GhostType &,const Array< UInt > &) const;
%ignore FEEngine::interpolateOnIntegrationPoints(const Array<Real> &,Array<Real> &,UInt,const ElementType&,const GhostType &) const;
-
+ %ignore FEEngine::onNodesAdded;
+ %ignore FEEngine::onNodesRemoved;
+ %ignore FEEngine::onElementsAdded;
+ %ignore FEEngine::onElementsChanged;
+ %ignore FEEngine::onElementsRemoved;
+ %ignore FEEngine::elementTypes;
}
%include "sparse_matrix.i"
%include "fe_engine.hh"
%rename(FreeBoundaryDirichlet) akantu::BC::Dirichlet::FreeBoundary;
%rename(FreeBoundaryNeumann) akantu::BC::Neumann::FreeBoundary;
%rename(PythonBoundary) akantu::BC::Dirichlet::PythonFunctor;
%include "boundary_condition_functor.hh"
%include "boundary_condition.hh"
%include "boundary_condition_python_functor.hh"
%include "communication_buffer.hh"
%include "data_accessor.hh"
-%include "synchronizer.hh"
-%include "synchronizer_registry.hh"
+//%include "synchronizer.hh"
+//%include "synchronizer_registry.hh"
%include "model.hh"
diff --git a/python/swig/solid_mechanics_model.i b/python/swig/solid_mechanics_model.i
index f9fe96049..492a77085 100644
--- a/python/swig/solid_mechanics_model.i
+++ b/python/swig/solid_mechanics_model.i
@@ -1,225 +1,204 @@
/**
* @file solid_mechanics_model.i
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Fabian Barras <fabian.barras@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Fri Dec 12 2014
* @date last modification: Wed Jan 06 2016
*
* @brief solid mechanics model wrapper
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
%{
#include "solid_mechanics_model.hh"
#include "sparse_matrix.hh"
#include "boundary_condition.hh"
#include "boundary_condition_functor.hh"
#include "boundary_condition_python_functor.hh"
#include "material_selector.hh"
#include "material_python.hh"
%}
namespace akantu {
%ignore SolidMechanicsModel::initFEEngineBoundary;
%ignore SolidMechanicsModel::initParallel;
%ignore SolidMechanicsModel::initArrays;
%ignore SolidMechanicsModel::initModel;
%ignore SolidMechanicsModel::initPBC;
ignore SolidMechanicsModel::initExplicit;
%ignore SolidMechanicsModel::isExplicit;
%ignore SolidMechanicsModel::updateCurrentPosition;
%ignore SolidMechanicsModel::updateAcceleration;
%ignore SolidMechanicsModel::updateIncrement;
%ignore SolidMechanicsModel::updatePreviousDisplacement;
%ignore SolidMechanicsModel::saveStressAndStrainBeforeDamage;
%ignore SolidMechanicsModel::updateEnergiesAfterDamage;
%ignore SolidMechanicsModel::solveLumped;
%ignore SolidMechanicsModel::explicitPred;
%ignore SolidMechanicsModel::explicitCorr;
%ignore SolidMechanicsModel::initSolver;
%ignore SolidMechanicsModel::initImplicit;
%ignore SolidMechanicsModel::initialAcceleration;
%ignore SolidMechanicsModel::testConvergence;
%ignore SolidMechanicsModel::testConvergenceIncrement;
%ignore SolidMechanicsModel::testConvergenceResidual;
%ignore SolidMechanicsModel::initVelocityDampingMatrix;
%ignore SolidMechanicsModel::getNbDataForElements;
%ignore SolidMechanicsModel::packElementData;
%ignore SolidMechanicsModel::unpackElementData;
%ignore SolidMechanicsModel::getNbDataToPack;
%ignore SolidMechanicsModel::getNbDataToUnpack;
%ignore SolidMechanicsModel::packData;
%ignore SolidMechanicsModel::unpackData;
%ignore SolidMechanicsModel::setMaterialSelector;
%ignore SolidMechanicsModel::getSolver;
%ignore SolidMechanicsModel::getSynchronizer;
%ignore Dumpable::registerExternalDumper;
%ignore Material::onNodesAdded;
%ignore Material::onNodesRemoved;
%ignore Material::onElementsAdded;
%ignore Material::onElementsRemoved;
%ignore Material::onElementsChanged;
}
%template(SolidMechanicsBoundaryCondition) akantu::BoundaryCondition<akantu::SolidMechanicsModel>;
%include "dumpable.hh"
print_self(SolidMechanicsModel)
%include "material.i"
%include "solid_mechanics_model.hh"
+%inline %{
+ namespace akantu{
+ void registerNewPythonMaterial(PyObject * obj, const akantu::ID & mat_type) {
+
+
+ MaterialFactory::getInstance().registerAllocator(
+ mat_type,
+ [obj](UInt, const ID &, SolidMechanicsModel & model,
+ const ID & id) -> std::unique_ptr<Material>
+ {
+ return std::make_unique<MaterialPython>(model, obj, id);
+ }
+ );
+ }
+ }
+%}
+
%extend akantu::SolidMechanicsModel {
/* ------------------------------------------------------------------------ */
void setPhysicalNamesMaterialSelector(){
akantu::MeshDataMaterialSelector<std::string> * selector = new
akantu::MeshDataMaterialSelector<std::string>("physical_names", *self);
self->setMaterialSelector(*selector);
}
/* ------------------------------------------------------------------------ */
- bool testConvergenceSccRes(Real tolerance) {
- Real error = 0;
- bool res = self->testConvergence<akantu::_scc_residual>(tolerance, error);
- return res;
+ void getResidual() {
+ AKANTU_EXCEPTION("Deprecated function. You should replace:\n"
+ "model.getResidual()\n"
+ "with\n"
+ "model.getInternalForce()\n");
}
- /* ------------------------------------------------------------------------ */
- void solveStaticDisplacement(Real tolerance, UInt max_iteration) {
- $self->solveStatic<akantu::_scm_newton_raphson_tangent_not_computed,
- akantu::_scc_residual>(tolerance, max_iteration);
- }
-
- /* ------------------------------------------------------------------------ */
- /// register an empty material of a given type
void registerNewPythonMaterial(PyObject * obj, const akantu::ID & mat_type) {
- std::pair<akantu::Parser::const_section_iterator,
- akantu::Parser::const_section_iterator>
- sub_sect = akantu::getStaticParser().getSubSections(akantu::_st_material);
-
- akantu::Parser::const_section_iterator it = sub_sect.first;
- for (; it != sub_sect.second; ++it) {
- if (it->getName() == mat_type) {
-
- AKANTU_DEBUG_ASSERT($self->materials_names_to_id.find(mat_type) ==
- $self->materials_names_to_id.end(),
- "A material with this name '"
- << mat_type << "' has already been registered. "
- << "Please use unique names for materials");
-
- UInt mat_count = $self->materials.size();
- $self->materials_names_to_id[mat_type] = mat_count;
-
- std::stringstream sstr_mat;
- sstr_mat << $self->getID() << ":" << mat_count << ":" << mat_type;
- akantu::ID mat_id = sstr_mat.str();
-
- akantu::Material * material = new akantu::MaterialPython(*$self, obj, mat_id);
- $self->materials.push_back(material);
-
- material->parseSection(*it);
- }
- }
+ AKANTU_EXCEPTION("Deprecated function. You should replace:\n"
+ "model.registerNewPythonMaterial(obj, mat_id)\n"
+ "with\n"
+ "akantu.registerNewPythonMaterial(obj, mat_id)\n\n"
+ "This MUST be done before instanciating the model.");
+
+ /* std::pair<akantu::Parser::const_section_iterator, */
+ /* akantu::Parser::const_section_iterator> */
+ /* sub_sect = akantu::getStaticParser().getSubSections(akantu::_st_material); */
+
+ /* akantu::Parser::const_section_iterator it = sub_sect.first; */
+ /* for (; it != sub_sect.second; ++it) { */
+ /* if (it->getName() == mat_type) { */
+
+ /* AKANTU_DEBUG_ASSERT($self->materials_names_to_id.find(mat_type) == */
+ /* $self->materials_names_to_id.end(), */
+ /* "A material with this name '" */
+ /* << mat_type << "' has already been registered. " */
+ /* << "Please use unique names for materials"); */
+
+ /* UInt mat_count = $self->materials.size(); */
+ /* $self->materials_names_to_id[mat_type] = mat_count; */
+
+ /* std::stringstream sstr_mat; */
+ /* sstr_mat << $self->getID() << ":" << mat_count << ":" << mat_type; */
+ /* akantu::ID mat_id = sstr_mat.str(); */
+
+ /* akantu::Material * material = new akantu::MaterialPython(*$self, obj, mat_id); */
+ /* $self->materials.push_back(material); */
+
+ /* material->parseSection(*it); */
+ /* } */
+ /* } */
}
/* ------------------------------------------------------------------------ */
void applyDirichletBC(PyObject * func_obj, const std::string & group_name) {
akantu::BC::PythonFunctorDirichlet functor(func_obj);
$self->applyBC(functor, group_name);
}
/* ------------------------------------------------------------------------ */
void applyNeumannBC(PyObject * func_obj, const std::string & group_name) {
akantu::BC::PythonFunctorNeumann functor(func_obj);
$self->applyBC(functor, group_name);
}
- /* ------------------------------------------------------------------------ */
- void solveDisplCorr(bool need_factorize, bool has_profile_changed) {
- akantu::Array<akantu::Real> & increment = $self->getIncrement();
-
- $self->solve<akantu::IntegrationScheme2ndOrder::_displacement_corrector>(
- increment, 1., need_factorize, has_profile_changed);
- }
-
- /* ------------------------------------------------------------------------ */
- void clearDispl() {
- akantu::Array<akantu::Real> & displ = $self->getDisplacement();
- displ.clear();
- }
-
- /* ------------------------------------------------------------------------ */
- void solveStep_TgModifIncr(Real tolerance, UInt max_iteration) {
- $self->solveStep<akantu::_scm_newton_raphson_tangent_modified,
- akantu::_scc_increment>(tolerance, max_iteration);
- }
-
- /* ------------------------------------------------------------------------ */
- void solveStep_TgIncr(Real tolerance, Real & error, UInt max_iteration) {
-
- $self->solveStep<akantu::_scm_newton_raphson_tangent,
- akantu::_scc_increment>(tolerance, error, max_iteration);
- }
-
- /* ------------------------------------------------------------------------ */
- void clearDisplVeloAcc() {
- akantu::Array<akantu::Real> & displ = $self->getDisplacement();
- akantu::Array<akantu::Real> & velo = $self->getVelocity();
- akantu::Array<akantu::Real> & acc = $self->getAcceleration();
-
- displ.clear();
- velo.clear();
- acc.clear();
- }
-
/* ------------------------------------------------------------------------ */
void applyUniformPressure(Real pressure, const std::string surface_name) {
UInt spatial_dimension = $self->getSpatialDimension();
akantu::Matrix<akantu::Real> surface_stress(spatial_dimension,
spatial_dimension, 0.0);
for (UInt i = 0; i < spatial_dimension; ++i) {
surface_stress(i, i) = -pressure;
}
$self->applyBC(akantu::BC::Neumann::FromStress(surface_stress),
surface_name);
}
/* ------------------------------------------------------------------------ */
void blockDOF(const std::string surface_name, SpacialDirection direction) {
$self->applyBC(akantu::BC::Dirichlet::FixedValue(0.0, direction),
surface_name);
}
}
+
diff --git a/python/swig/sparse_matrix.i b/python/swig/sparse_matrix.i
index 50acfc688..f3c2dc6b0 100644
--- a/python/swig/sparse_matrix.i
+++ b/python/swig/sparse_matrix.i
@@ -1,34 +1,34 @@
%include "sparse_matrix.hh"
%pythoncode %{
import scipy.sparse
import numpy as _np
class AkantuSparseMatrix (scipy.sparse.coo_matrix) :
def __init__(self,aka_sparse):
self.aka_sparse = aka_sparse
matrix_type = self.aka_sparse.getSparseMatrixType()
- sz = self.aka_sparse.getSize()
+ sz = self.aka_sparse.size()
row = self.aka_sparse.getIRN()[:,0] -1
col = self.aka_sparse.getJCN()[:,0] -1
data = self.aka_sparse.getA()[:,0]
row = row.copy()
col = col.copy()
data = data.copy()
if matrix_type == _symmetric:
non_diags = (row != col)
row_sup = col[non_diags]
col_sup = row[non_diags]
data_sup = data[non_diags]
col = _np.concatenate((col,col_sup))
row = _np.concatenate((row,row_sup))
data = _np.concatenate((data,data_sup))
scipy.sparse.coo_matrix.__init__(self,(data, (row,col)),shape=(sz,sz))
%}
diff --git a/src/common/aka_array_tmpl.hh b/src/common/aka_array_tmpl.hh
index 90c468ab8..97fefce38 100644
--- a/src/common/aka_array_tmpl.hh
+++ b/src/common/aka_array_tmpl.hh
@@ -1,1275 +1,1279 @@
/**
* @file aka_array_tmpl.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Thu Jul 15 2010
* @date last modification: Fri Jan 22 2016
*
* @brief Inline functions of the classes Array<T> and ArrayBase
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
/* Inline Functions Array<T> */
/* -------------------------------------------------------------------------- */
#include "aka_array.hh"
/* -------------------------------------------------------------------------- */
#include <memory>
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_AKA_ARRAY_TMPL_HH__
#define __AKANTU_AKA_ARRAY_TMPL_HH__
namespace akantu {
namespace debug {
struct ArrayException : public Exception {};
} // namespace debug
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
inline T & Array<T, is_scal>::operator()(UInt i, UInt j) {
AKANTU_DEBUG_ASSERT(size_ > 0, "The array \"" << id << "\" is empty");
AKANTU_DEBUG_ASSERT((i < size_) && (j < nb_component),
"The value at position ["
<< i << "," << j << "] is out of range in array \""
<< id << "\"");
return values[i * nb_component + j];
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
inline const T & Array<T, is_scal>::operator()(UInt i, UInt j) const {
AKANTU_DEBUG_ASSERT(size_ > 0, "The array \"" << id << "\" is empty");
AKANTU_DEBUG_ASSERT((i < size_) && (j < nb_component),
"The value at position ["
<< i << "," << j << "] is out of range in array \""
<< id << "\"");
return values[i * nb_component + j];
}
template <class T, bool is_scal>
inline T & Array<T, is_scal>::operator[](UInt i) {
AKANTU_DEBUG_ASSERT(size_ > 0, "The array \"" << id << "\" is empty");
AKANTU_DEBUG_ASSERT((i < size_ * nb_component),
"The value at position ["
<< i << "] is out of range in array \"" << id
<< "\"");
return values[i];
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
inline const T & Array<T, is_scal>::operator[](UInt i) const {
AKANTU_DEBUG_ASSERT(size_ > 0, "The array \"" << id << "\" is empty");
AKANTU_DEBUG_ASSERT((i < size_ * nb_component),
"The value at position ["
<< i << "] is out of range in array \"" << id
<< "\"");
return values[i];
}
/* -------------------------------------------------------------------------- */
/**
* append a tuple to the array with the value value for all components
* @param value the new last tuple or the array will contain nb_component copies
* of value
*/
template <class T, bool is_scal>
inline void Array<T, is_scal>::push_back(const T & value) {
resizeUnitialized(size_ + 1, true, value);
}
/* -------------------------------------------------------------------------- */
/**
* append a tuple to the array
* @param new_elem a C-array containing the values to be copied to the end of
* the array */
// template <class T, bool is_scal>
// inline void Array<T, is_scal>::push_back(const T new_elem[]) {
// UInt pos = size_;
// resizeUnitialized(size_ + 1, false);
// T * tmp = values + nb_component * pos;
// std::uninitialized_copy(new_elem, new_elem + nb_component, tmp);
// }
/* -------------------------------------------------------------------------- */
#ifndef SWIG
/**
* append a matrix or a vector to the array
* @param new_elem a reference to a Matrix<T> or Vector<T> */
template <class T, bool is_scal>
template <template <typename> class C, typename>
inline void Array<T, is_scal>::push_back(const C<T> & new_elem) {
AKANTU_DEBUG_ASSERT(
nb_component == new_elem.size(),
"The vector("
<< new_elem.size()
<< ") as not a size compatible with the Array (nb_component="
<< nb_component << ").");
UInt pos = size_;
resizeUnitialized(size_ + 1, false);
T * tmp = values + nb_component * pos;
std::uninitialized_copy(new_elem.storage(), new_elem.storage() + nb_component,
tmp);
}
/* -------------------------------------------------------------------------- */
/**
* append a tuple to the array
* @param it an iterator to the tuple to be copied to the end of the array */
template <class T, bool is_scal>
template <class Ret>
inline void
Array<T, is_scal>::push_back(const Array<T, is_scal>::iterator<Ret> & it) {
UInt pos = size_;
resizeUnitialized(size_ + 1, false);
T * tmp = values + nb_component * pos;
T * new_elem = it.data();
std::uninitialized_copy(new_elem, new_elem + nb_component, tmp);
}
#endif
/* -------------------------------------------------------------------------- */
/**
* erase an element. If the erased element is not the last of the array, the
* last element is moved into the hole in order to maintain contiguity. This
* may invalidate existing iterators (For instance an iterator obtained by
* Array::end() is no longer correct) and will change the order of the
* elements.
* @param i index of element to erase
*/
template <class T, bool is_scal> inline void Array<T, is_scal>::erase(UInt i) {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT((size_ > 0), "The array is empty");
AKANTU_DEBUG_ASSERT((i < size_),
"The element at position [" << i << "] is out of range ("
<< i << ">=" << size_ << ")");
if (i != (size_ - 1)) {
for (UInt j = 0; j < nb_component; ++j) {
values[i * nb_component + j] = values[(size_ - 1) * nb_component + j];
}
}
resize(size_ - 1);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Subtract another array entry by entry from this array in place. Both arrays
* must
* have the same size and nb_component. If the arrays have different shapes,
* code compiled in debug mode will throw an expeption and optimised code
* will behave in an unpredicted manner
* @param other array to subtract from this
* @return reference to modified this
*/
template <class T, bool is_scal>
Array<T, is_scal> & Array<T, is_scal>::
operator-=(const Array<T, is_scal> & vect) {
AKANTU_DEBUG_ASSERT((size_ == vect.size_) &&
(nb_component == vect.nb_component),
"The too array don't have the same sizes");
T * a = values;
T * b = vect.storage();
for (UInt i = 0; i < size_ * nb_component; ++i) {
*a -= *b;
++a;
++b;
}
return *this;
}
/* -------------------------------------------------------------------------- */
/**
* Add another array entry by entry to this array in place. Both arrays must
* have the same size and nb_component. If the arrays have different shapes,
* code compiled in debug mode will throw an expeption and optimised code
* will behave in an unpredicted manner
* @param other array to add to this
* @return reference to modified this
*/
template <class T, bool is_scal>
Array<T, is_scal> & Array<T, is_scal>::
operator+=(const Array<T, is_scal> & vect) {
- AKANTU_DEBUG_ASSERT((size_ == vect.size) &&
+ AKANTU_DEBUG_ASSERT((size_ == vect.size()) &&
(nb_component == vect.nb_component),
"The too array don't have the same sizes");
T * a = values;
T * b = vect.storage();
for (UInt i = 0; i < size_ * nb_component; ++i) {
*a++ += *b++;
}
return *this;
}
/* -------------------------------------------------------------------------- */
/**
* Multiply all entries of this array by a scalar in place
* @param alpha scalar multiplicant
* @return reference to modified this
*/
+
+#ifndef SWIG
template <class T, bool is_scal>
Array<T, is_scal> & Array<T, is_scal>::operator*=(const T & alpha) {
T * a = values;
for (UInt i = 0; i < size_ * nb_component; ++i) {
*a++ *= alpha;
}
return *this;
}
+#endif
+
/* -------------------------------------------------------------------------- */
/**
* Compare this array element by element to another.
* @param other array to compare to
* @return true it all element are equal and arrays have the same shape, else
* false
*/
template <class T, bool is_scal>
bool Array<T, is_scal>::operator==(const Array<T, is_scal> & array) const {
bool equal = nb_component == array.nb_component && size_ == array.size_ &&
id == array.id;
if (!equal)
return false;
if (values == array.storage())
return true;
else
return std::equal(values, values + size_ * nb_component, array.storage());
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
bool Array<T, is_scal>::operator!=(const Array<T, is_scal> & array) const {
return !operator==(array);
}
/* -------------------------------------------------------------------------- */
#ifndef SWIG
/**
* set all tuples of the array to a given vector or matrix
* @param vm Matrix or Vector to fill the array with
*/
template <class T, bool is_scal>
template <template <typename> class C, typename>
inline void Array<T, is_scal>::set(const C<T> & vm) {
AKANTU_DEBUG_ASSERT(
nb_component == vm.size(),
"The size of the object does not match the number of components");
for (T * it = values; it < values + nb_component * size_;
it += nb_component) {
std::copy(vm.storage(), vm.storage() + nb_component, it);
}
}
#endif
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
void Array<T, is_scal>::append(const Array<T> & other) {
AKANTU_DEBUG_ASSERT(
nb_component == other.nb_component,
"Cannot append an array with a different number of component");
UInt old_size = this->size_;
this->resizeUnitialized(this->size_ + other.size(), false);
T * tmp = values + nb_component * old_size;
std::uninitialized_copy(other.storage(),
other.storage() + other.size() * nb_component, tmp);
}
/* -------------------------------------------------------------------------- */
/* Functions Array<T, is_scal> */
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
Array<T, is_scal>::Array(UInt size, UInt nb_component, const ID & id)
: ArrayBase(id), values(nullptr) {
AKANTU_DEBUG_IN();
allocate(size, nb_component);
if (!is_scal) {
T val = T();
std::uninitialized_fill(values, values + size * nb_component, val);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
Array<T, is_scal>::Array(UInt size, UInt nb_component, const T def_values[],
const ID & id)
: ArrayBase(id), values(NULL) {
AKANTU_DEBUG_IN();
allocate(size, nb_component);
T * tmp = values;
for (UInt i = 0; i < size; ++i) {
tmp = values + nb_component * i;
std::uninitialized_copy(def_values, def_values + nb_component, tmp);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
Array<T, is_scal>::Array(UInt size, UInt nb_component, const T & value,
const ID & id)
: ArrayBase(id), values(nullptr) {
AKANTU_DEBUG_IN();
allocate(size, nb_component);
std::uninitialized_fill_n(values, size * nb_component, value);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
Array<T, is_scal>::Array(const Array<T, is_scal> & vect, bool deep,
const ID & id)
: ArrayBase(vect) {
AKANTU_DEBUG_IN();
this->id = (id == "") ? vect.id : id;
if (deep) {
allocate(vect.size_, vect.nb_component);
T * tmp = values;
std::uninitialized_copy(vect.storage(),
vect.storage() + size_ * nb_component, tmp);
} else {
this->values = vect.storage();
this->size_ = vect.size_;
this->nb_component = vect.nb_component;
this->allocated_size = vect.allocated_size;
this->size_of_type = vect.size_of_type;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
#ifndef SWIG
template <class T, bool is_scal>
Array<T, is_scal>::Array(const std::vector<T> & vect) {
AKANTU_DEBUG_IN();
this->id = "";
allocate(vect.size(), 1);
T * tmp = values;
std::uninitialized_copy(&(vect[0]), &(vect[size_ - 1]), tmp);
AKANTU_DEBUG_OUT();
}
#endif
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal> Array<T, is_scal>::~Array() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG(dblAccessory,
"Freeing " << printMemorySize<T>(allocated_size * nb_component)
<< " (" << id << ")");
if (values) {
if (!is_scal)
for (UInt i = 0; i < size_ * nb_component; ++i) {
T * obj = values + i;
obj->~T();
}
free(values);
}
size_ = allocated_size = 0;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* perform the allocation for the constructors
* @param size is the size of the array
* @param nb_component is the number of component of the array
*/
template <class T, bool is_scal>
void Array<T, is_scal>::allocate(UInt size, UInt nb_component) {
AKANTU_DEBUG_IN();
if (size == 0) {
values = nullptr;
} else {
values = static_cast<T *>(malloc(nb_component * size * sizeof(T)));
AKANTU_DEBUG_ASSERT(values != nullptr,
"Cannot allocate "
<< printMemorySize<T>(size * nb_component) << " ("
<< id << ")");
}
if (values == nullptr) {
this->size_ = this->allocated_size = 0;
} else {
AKANTU_DEBUG(dblAccessory,
"Allocated " << printMemorySize<T>(size * nb_component) << " ("
<< id << ")");
this->size_ = this->allocated_size = size;
}
this->size_of_type = sizeof(T);
this->nb_component = nb_component;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
void Array<T, is_scal>::reserve(UInt new_size) {
UInt tmp_size = this->size_;
resizeUnitialized(new_size, false);
this->size_ = tmp_size;
}
/* -------------------------------------------------------------------------- */
/**
* change the size of the array and allocate or free memory if needed. If the
* size increases, the new tuples are filled with zeros
* @param new_size new number of tuples contained in the array */
template <class T, bool is_scal> void Array<T, is_scal>::resize(UInt new_size) {
resizeUnitialized(new_size, !is_scal);
}
/* -------------------------------------------------------------------------- */
/**
* change the size of the array and allocate or free memory if needed. If the
* size increases, the new tuples are filled with zeros
* @param new_size new number of tuples contained in the array */
template <class T, bool is_scal>
void Array<T, is_scal>::resize(UInt new_size, const T & val) {
this->resizeUnitialized(new_size, true, val);
}
/* -------------------------------------------------------------------------- */
/**
* change the size of the array and allocate or free memory if needed.
* @param new_size new number of tuples contained in the array */
template <class T, bool is_scal>
void Array<T, is_scal>::resizeUnitialized(UInt new_size, bool fill,
const T & val) {
// AKANTU_DEBUG_IN();
// free some memory
if (new_size <= allocated_size) {
if (!is_scal) {
T * old_values = values;
if (new_size < size_) {
for (UInt i = new_size * nb_component; i < size_ * nb_component; ++i) {
T * obj = old_values + i;
obj->~T();
}
}
}
if (allocated_size - new_size > AKANTU_MIN_ALLOCATION) {
AKANTU_DEBUG(dblAccessory,
"Freeing " << printMemorySize<T>((allocated_size - size_) *
nb_component)
<< " (" << id << ")");
// Normally there are no allocation problem when reducing an array
if (new_size == 0) {
free(values);
values = nullptr;
} else {
auto * tmp_ptr = static_cast<T *>(
realloc(values, new_size * nb_component * sizeof(T)));
if (tmp_ptr == nullptr) {
AKANTU_EXCEPTION("Cannot free data ("
<< id << ")"
<< " [current allocated size : " << allocated_size
<< " | "
<< "requested size : " << new_size << "]");
}
values = tmp_ptr;
}
allocated_size = new_size;
}
} else {
// allocate more memory
UInt size_to_alloc = (new_size - allocated_size < AKANTU_MIN_ALLOCATION)
? allocated_size + AKANTU_MIN_ALLOCATION
: new_size;
auto * tmp_ptr = static_cast<T *>(
realloc(values, size_to_alloc * nb_component * sizeof(T)));
AKANTU_DEBUG_ASSERT(
tmp_ptr != nullptr,
"Cannot allocate " << printMemorySize<T>(size_to_alloc * nb_component));
if (tmp_ptr == nullptr) {
AKANTU_DEBUG_ERROR("Cannot allocate more data ("
<< id << ")"
<< " [current allocated size : " << allocated_size
<< " | "
<< "requested size : " << new_size << "]");
}
AKANTU_DEBUG(dblAccessory,
"Allocating " << printMemorySize<T>(
(size_to_alloc - allocated_size) * nb_component));
allocated_size = size_to_alloc;
values = tmp_ptr;
}
if (fill && this->size_ < new_size) {
std::uninitialized_fill(values + size_ * nb_component,
values + new_size * nb_component, val);
}
size_ = new_size;
// AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* change the number of components by interlacing data
* @param multiplicator number of interlaced components add
* @param block_size blocks of data in the array
* Examaple for block_size = 2, multiplicator = 2
* array = oo oo oo -> new array = oo nn nn oo nn nn oo nn nn */
template <class T, bool is_scal>
void Array<T, is_scal>::extendComponentsInterlaced(UInt multiplicator,
UInt block_size) {
AKANTU_DEBUG_IN();
if (multiplicator == 1)
return;
AKANTU_DEBUG_ASSERT(multiplicator > 1, "invalid multiplicator");
AKANTU_DEBUG_ASSERT(nb_component % block_size == 0,
"stride must divide actual number of components");
values = static_cast<T *>(
realloc(values, nb_component * multiplicator * size_ * sizeof(T)));
UInt new_component = nb_component / block_size * multiplicator;
for (UInt i = 0, k = size_ - 1; i < size_; ++i, --k) {
for (UInt j = 0; j < new_component; ++j) {
UInt m = new_component - j - 1;
UInt n = m / multiplicator;
for (UInt l = 0, p = block_size - 1; l < block_size; ++l, --p) {
values[k * nb_component * multiplicator + m * block_size + p] =
values[k * nb_component + n * block_size + p];
}
}
}
nb_component = nb_component * multiplicator;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* search elem in the array, return the position of the first occurrence or
* -1 if not found
* @param elem the element to look for
* @return index of the first occurrence of elem or -1 if elem is not present
*/
template <class T, bool is_scal>
UInt Array<T, is_scal>::find(const T & elem) const {
AKANTU_DEBUG_IN();
auto begin = this->begin();
auto end = this->end();
auto it = std::find(begin, end, elem);
AKANTU_DEBUG_OUT();
return (it != end) ? it - begin : UInt(-1);
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal> UInt Array<T, is_scal>::find(T elem[]) const {
AKANTU_DEBUG_IN();
T * it = values;
UInt i = 0;
for (; i < size_; ++i) {
if (*it == elem[0]) {
T * cit = it;
UInt c = 0;
for (; (c < nb_component) && (*cit == elem[c]); ++c, ++cit)
;
if (c == nb_component) {
AKANTU_DEBUG_OUT();
return i;
}
}
it += nb_component;
}
return UInt(-1);
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
template <template <typename> class C, typename>
inline UInt Array<T, is_scal>::find(const C<T> & elem) {
AKANTU_DEBUG_ASSERT(elem.size() == nb_component,
"Cannot find an element with a wrong size ("
<< elem.size() << ") != " << nb_component);
return this->find(elem.storage());
}
/* -------------------------------------------------------------------------- */
/**
* copy the content of another array. This overwrites the current content.
* @param other Array to copy into this array. It has to have the same
* nb_component as this. If compiled in debug mode, an incorrect other will
* result in an exception being thrown. Optimised code may result in
* unpredicted behaviour.
*/
template <class T, bool is_scal>
void Array<T, is_scal>::copy(const Array<T, is_scal> & vect,
bool no_sanity_check) {
AKANTU_DEBUG_IN();
if (!no_sanity_check)
if (vect.nb_component != nb_component)
AKANTU_DEBUG_ERROR(
"The two arrays do not have the same number of components");
resize((vect.size_ * vect.nb_component) / nb_component);
T * tmp = values;
std::uninitialized_copy(vect.storage(), vect.storage() + size_ * nb_component,
tmp);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <bool is_scal> class ArrayPrintHelper {
public:
template <typename T>
static void print_content(const Array<T> & vect, std::ostream & stream,
int indent) {
if (AKANTU_DEBUG_TEST(dblDump) || AKANTU_DEBUG_LEVEL_IS_TEST()) {
std::string space;
for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
;
stream << space << " + values : {";
for (UInt i = 0; i < vect.size(); ++i) {
stream << "{";
for (UInt j = 0; j < vect.getNbComponent(); ++j) {
stream << vect(i, j);
if (j != vect.getNbComponent() - 1)
stream << ", ";
}
stream << "}";
if (i != vect.size() - 1)
stream << ", ";
}
stream << "}" << std::endl;
}
}
};
template <> class ArrayPrintHelper<false> {
public:
template <typename T>
static void print_content(__attribute__((unused)) const Array<T> & vect,
__attribute__((unused)) std::ostream & stream,
__attribute__((unused)) int indent) {}
};
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
void Array<T, is_scal>::printself(std::ostream & stream, int indent) const {
std::string space;
for (Int i = 0; i < indent; i++, space += AKANTU_INDENT)
;
std::streamsize prec = stream.precision();
std::ios_base::fmtflags ff = stream.flags();
stream.setf(std::ios_base::showbase);
stream.precision(2);
stream << space << "Array<" << debug::demangle(typeid(T).name()) << "> ["
<< std::endl;
stream << space << " + id : " << this->id << std::endl;
stream << space << " + size : " << this->size_ << std::endl;
stream << space << " + nb_component : " << this->nb_component << std::endl;
stream << space << " + allocated size : " << this->allocated_size
<< std::endl;
stream << space << " + memory size : "
<< printMemorySize<T>(allocated_size * nb_component) << std::endl;
if (!AKANTU_DEBUG_LEVEL_IS_TEST())
stream << space << " + address : " << std::hex << this->values
<< std::dec << std::endl;
stream.precision(prec);
stream.flags(ff);
ArrayPrintHelper<is_scal>::print_content(*this, stream, indent);
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
/* Inline Functions ArrayBase */
/* -------------------------------------------------------------------------- */
inline UInt ArrayBase::getMemorySize() const {
return allocated_size * nb_component * size_of_type;
}
inline void ArrayBase::empty() { size_ = 0; }
#ifndef SWIG
/* -------------------------------------------------------------------------- */
/* Iterators */
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
template <class R, class daughter, class IR, bool is_tensor>
class Array<T, is_scal>::iterator_internal {
public:
using value_type = R;
using pointer = R *;
using reference = R &;
using proxy = typename R::proxy;
using const_proxy = const typename R::proxy;
using const_reference = const R &;
using internal_value_type = IR;
using internal_pointer = IR *;
using difference_type = std::ptrdiff_t;
using iterator_category = std::random_access_iterator_tag;
public:
iterator_internal() = default;
iterator_internal(pointer_type data, UInt _offset)
: _offset(_offset), initial(data), ret(nullptr), ret_ptr(data) {
AKANTU_DEBUG_ERROR(
"The constructor should never be called it is just an ugly trick...");
}
iterator_internal(std::unique_ptr<internal_value_type> && wrapped)
: _offset(wrapped->size()), initial(wrapped->storage()),
ret(std::move(wrapped)), ret_ptr(ret->storage()) {}
iterator_internal(const iterator_internal & it) {
if (this != &it) {
this->_offset = it._offset;
this->initial = it.initial;
this->ret_ptr = it.ret_ptr;
this->ret = std::make_unique<internal_value_type>(*it.ret, false);
}
}
iterator_internal(iterator_internal && it) = default;
virtual ~iterator_internal() = default;
inline iterator_internal & operator=(const iterator_internal & it) {
if (this != &it) {
this->_offset = it._offset;
this->initial = it.initial;
this->ret_ptr = it.ret_ptr;
if (this->ret)
this->ret->shallowCopy(*it.ret);
else
this->ret = std::make_unique<internal_value_type>(*it.ret, false);
}
return *this;
}
UInt getCurrentIndex() {
return (this->ret_ptr - this->initial) / this->_offset;
};
inline reference operator*() {
ret->values = ret_ptr;
return *ret;
};
inline const_reference operator*() const {
ret->values = ret_ptr;
return *ret;
};
inline pointer operator->() {
ret->values = ret_ptr;
return ret.get();
};
inline daughter & operator++() {
ret_ptr += _offset;
return static_cast<daughter &>(*this);
};
inline daughter & operator--() {
ret_ptr -= _offset;
return static_cast<daughter &>(*this);
};
inline daughter & operator+=(const UInt n) {
ret_ptr += _offset * n;
return static_cast<daughter &>(*this);
}
inline daughter & operator-=(const UInt n) {
ret_ptr -= _offset * n;
return static_cast<daughter &>(*this);
}
inline proxy operator[](const UInt n) {
ret->values = ret_ptr + n * _offset;
return proxy(*ret);
}
inline const_proxy operator[](const UInt n) const {
ret->values = ret_ptr + n * _offset;
return const_proxy(*ret);
}
inline bool operator==(const iterator_internal & other) const {
return this->ret_ptr == other.ret_ptr;
}
inline bool operator!=(const iterator_internal & other) const {
return this->ret_ptr != other.ret_ptr;
}
inline bool operator<(const iterator_internal & other) const {
return this->ret_ptr < other.ret_ptr;
}
inline bool operator<=(const iterator_internal & other) const {
return this->ret_ptr <= other.ret_ptr;
}
inline bool operator>(const iterator_internal & other) const {
return this->ret_ptr > other.ret_ptr;
}
inline bool operator>=(const iterator_internal & other) const {
return this->ret_ptr >= other.ret_ptr;
}
inline daughter operator+(difference_type n) {
daughter tmp(static_cast<daughter &>(*this));
tmp += n;
return tmp;
}
inline daughter operator-(difference_type n) {
daughter tmp(static_cast<daughter &>(*this));
tmp -= n;
return tmp;
}
inline difference_type operator-(const iterator_internal & b) {
return (this->ret_ptr - b.ret_ptr) / _offset;
}
inline pointer_type data() const { return ret_ptr; }
inline difference_type offset() const { return _offset; }
protected:
UInt _offset{0};
pointer_type initial{nullptr};
std::unique_ptr<internal_value_type> ret{nullptr};
pointer_type ret_ptr{nullptr};
};
/* -------------------------------------------------------------------------- */
/**
* Specialization for scalar types
*/
template <class T, bool is_scal>
template <class R, class daughter, class IR>
class Array<T, is_scal>::iterator_internal<R, daughter, IR, false> {
public:
using value_type = R;
using pointer = R *;
using reference = R &;
using const_reference = const R &;
using internal_value_type = IR;
using internal_pointer = IR *;
using difference_type = std::ptrdiff_t;
using iterator_category = std::random_access_iterator_tag;
public:
iterator_internal(pointer data = nullptr) : ret(data), initial(data){};
iterator_internal(const iterator_internal & it) = default;
iterator_internal(iterator_internal && it) = default;
virtual ~iterator_internal() = default;
inline iterator_internal & operator=(const iterator_internal & it) = default;
UInt getCurrentIndex() { return (this->ret - this->initial); };
inline reference operator*() { return *ret; };
inline const_reference operator*() const { return *ret; };
inline pointer operator->() { return ret; };
inline daughter & operator++() {
++ret;
return static_cast<daughter &>(*this);
};
inline daughter & operator--() {
--ret;
return static_cast<daughter &>(*this);
};
inline daughter & operator+=(const UInt n) {
ret += n;
return static_cast<daughter &>(*this);
}
inline daughter & operator-=(const UInt n) {
ret -= n;
return static_cast<daughter &>(*this);
}
inline reference operator[](const UInt n) { return ret[n]; }
inline bool operator==(const iterator_internal & other) const {
return ret == other.ret;
}
inline bool operator!=(const iterator_internal & other) const {
return ret != other.ret;
}
inline bool operator<(const iterator_internal & other) const {
return ret < other.ret;
}
inline bool operator<=(const iterator_internal & other) const {
return ret <= other.ret;
}
inline bool operator>(const iterator_internal & other) const {
return ret > other.ret;
}
inline bool operator>=(const iterator_internal & other) const {
return ret >= other.ret;
}
inline daughter operator-(difference_type n) { return daughter(ret - n); }
inline daughter operator+(difference_type n) { return daughter(ret + n); }
inline difference_type operator-(const iterator_internal & b) {
return ret - b.ret;
}
inline pointer data() const { return ret; }
protected:
pointer ret{nullptr};
pointer initial{nullptr};
};
/* -------------------------------------------------------------------------- */
/* Begin/End functions implementation */
/* -------------------------------------------------------------------------- */
namespace detail {
template <class Tuple, size_t... Is>
constexpr auto take_front_impl(Tuple && t, std::index_sequence<Is...>) {
return std::make_tuple(std::get<Is>(std::forward<Tuple>(t))...);
}
template <size_t N, class Tuple> constexpr auto take_front(Tuple && t) {
return take_front_impl(std::forward<Tuple>(t),
std::make_index_sequence<N>{});
}
template <typename... V> constexpr auto product_all(V &&... v) {
std::common_type_t<int, V...> result = 1;
(void)std::initializer_list<int>{(result *= v, 0)...};
return result;
}
template <typename... T> std::string to_string_all(T &&... t) {
if (sizeof...(T) == 0)
return "";
std::stringstream ss;
bool noComma = true;
ss << "(";
(void)std::initializer_list<bool>{
(ss << (noComma ? "" : ", ") << t, noComma = false)...};
ss << ")";
return ss.str();
}
template <std::size_t N> struct InstantiationHelper {
template <typename type, typename T, typename... Ns>
static auto instantiate(T && data, Ns... ns) {
return std::make_unique<type>(data, ns...);
}
};
template <> struct InstantiationHelper<0> {
template <typename type, typename T> static auto instantiate(T && data) {
return data;
}
};
template <typename Arr, typename T, typename... Ns>
decltype(auto) get_iterator(Arr && array, T * data, Ns &&... ns) {
using type = IteratorHelper_t<sizeof...(Ns) - 1, T>;
using array_type = std::decay_t<Arr>;
using iterator =
std::conditional_t<std::is_const<std::remove_reference_t<Arr>>::value,
typename array_type::template const_iterator<type>,
typename array_type::template iterator<type>>;
static_assert(sizeof...(Ns), "You should provide a least one size");
if (array.getNbComponent() * array.size() !=
product_all(std::forward<Ns>(ns)...)) {
AKANTU_CUSTOM_EXCEPTION_INFO(
debug::ArrayException(),
"The iterator on "
<< debug::demangle(typeid(Arr).name())
<< to_string_all(array.size(), array.getNbComponent())
<< "is not compatible with the type "
<< debug::demangle(typeid(type).name()) << to_string_all(ns...));
}
auto && wrapped = aka::apply(
[&](auto... n) {
return InstantiationHelper<sizeof...(n)>::template instantiate<type>(
data, n...);
},
take_front<sizeof...(Ns) - 1>(std::make_tuple(ns...)));
return iterator(std::move(wrapped));
}
} // namespace detail
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
template <typename... Ns>
inline decltype(auto) Array<T, is_scal>::begin(Ns &&... ns) {
return detail::get_iterator(*this, values, std::forward<Ns>(ns)..., size_);
}
template <class T, bool is_scal>
template <typename... Ns>
inline decltype(auto) Array<T, is_scal>::end(Ns &&... ns) {
return detail::get_iterator(*this, values + nb_component * size_,
std::forward<Ns>(ns)..., size_);
}
template <class T, bool is_scal>
template <typename... Ns>
inline decltype(auto) Array<T, is_scal>::begin(Ns &&... ns) const {
return detail::get_iterator(*this, values, std::forward<Ns>(ns)..., size_);
}
template <class T, bool is_scal>
template <typename... Ns>
inline decltype(auto) Array<T, is_scal>::end(Ns &&... ns) const {
return detail::get_iterator(*this, values + nb_component * size_,
std::forward<Ns>(ns)..., size_);
}
template <class T, bool is_scal>
template <typename... Ns>
inline decltype(auto) Array<T, is_scal>::begin_reinterpret(Ns &&... ns) {
return detail::get_iterator(*this, values, std::forward<Ns>(ns)...);
}
template <class T, bool is_scal>
template <typename... Ns>
inline decltype(auto) Array<T, is_scal>::end_reinterpret(Ns &&... ns) {
return detail::get_iterator(
*this, values + detail::product_all(std::forward<Ns>(ns)...),
std::forward<Ns>(ns)...);
}
template <class T, bool is_scal>
template <typename... Ns>
inline decltype(auto) Array<T, is_scal>::begin_reinterpret(Ns &&... ns) const {
return detail::get_iterator(*this, values, std::forward<Ns>(ns)...);
}
template <class T, bool is_scal>
template <typename... Ns>
inline decltype(auto) Array<T, is_scal>::end_reinterpret(Ns &&... ns) const {
return detail::get_iterator(
*this, values + detail::product_all(std::forward<Ns>(ns)...),
std::forward<Ns>(ns)...);
}
/* -------------------------------------------------------------------------- */
/* Views */
/* -------------------------------------------------------------------------- */
namespace detail {
template <typename Array, typename... Ns> class ArrayView {
using tuple = std::tuple<Ns...>;
public:
ArrayView(Array && array, Ns... ns)
: array(std::forward<Array>(array)), sizes(std::move(ns)...) {}
ArrayView(ArrayView && array_view) = default;
ArrayView & operator=(const ArrayView & array_view) = default;
ArrayView & operator=(ArrayView && array_view) = default;
decltype(auto) begin() {
return aka::apply(
[&](auto &&... ns) { return array.begin_reinterpret(ns...); }, sizes);
}
decltype(auto) begin() const {
return aka::apply(
[&](auto &&... ns) { return array.begin_reinterpret(ns...); }, sizes);
}
decltype(auto) end() {
return aka::apply(
[&](auto &&... ns) { return array.end_reinterpret(ns...); }, sizes);
}
decltype(auto) end() const {
return aka::apply(
[&](auto &&... ns) { return array.end_reinterpret(ns...); }, sizes);
}
decltype(auto) size() const {
return std::get<std::tuple_size<tuple>::value - 1>(sizes);
}
decltype(auto) dims() const { return std::tuple_size<tuple>::value - 1; }
private:
Array array;
tuple sizes;
};
} // namespace detail
/* -------------------------------------------------------------------------- */
template <typename Array, typename... Ns>
decltype(auto) make_view(Array && array, Ns... ns) {
auto size = std::forward<decltype(array)>(array).size() *
std::forward<decltype(array)>(array).getNbComponent() /
detail::product_all(ns...);
return detail::ArrayView<Array, std::common_type_t<size_t, Ns>...,
std::common_type_t<size_t, decltype(size)>>(
std::forward<Array>(array), std::move(ns)..., size);
}
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
template <typename R>
class Array<T, is_scal>::const_iterator
: public iterator_internal<const R, Array<T, is_scal>::const_iterator<R>,
R> {
public:
using parent = iterator_internal<const R, const_iterator, R>;
using value_type = typename parent::value_type;
using pointer = typename parent::pointer;
using reference = typename parent::reference;
using difference_type = typename parent::difference_type;
using iterator_category = typename parent::iterator_category;
public:
const_iterator() : parent(){};
// const_iterator(pointer_type data, UInt offset) : parent(data, offset) {}
// const_iterator(pointer warped) : parent(warped) {}
// const_iterator(const parent & it) : parent(it) {}
const_iterator(const const_iterator & it) = default;
const_iterator(const_iterator && it) = default;
template <typename P, typename = std::enable_if_t<not is_tensor<P>::value>>
const_iterator(P * data) : parent(data) {}
template <typename UP_P,
typename =
std::enable_if_t<is_tensor<typename UP_P::element_type>::value>>
const_iterator(UP_P && tensor) : parent(std::forward<UP_P>(tensor)) {}
const_iterator & operator=(const const_iterator & it) = default;
};
template <class T, class R, bool is_tensor_ = is_tensor<R>::value>
struct ConstConverterIteratorHelper {
using const_iterator = typename Array<T>::template const_iterator<R>;
using iterator = typename Array<T>::template iterator<R>;
static inline const_iterator convert(const iterator & it) {
return const_iterator(std::unique_ptr<R>(new R(*it, false)));
}
};
template <class T, class R> struct ConstConverterIteratorHelper<T, R, false> {
using const_iterator = typename Array<T>::template const_iterator<R>;
using iterator = typename Array<T>::template iterator<R>;
static inline const_iterator convert(const iterator & it) {
return const_iterator(it.data());
}
};
template <class T, bool is_scal>
template <typename R>
class Array<T, is_scal>::iterator
: public iterator_internal<R, Array<T, is_scal>::iterator<R>> {
public:
using parent = iterator_internal<R, iterator>;
using value_type = typename parent::value_type;
using pointer = typename parent::pointer;
using reference = typename parent::reference;
using difference_type = typename parent::difference_type;
using iterator_category = typename parent::iterator_category;
public:
iterator() : parent(){};
iterator(const iterator & it) = default;
iterator(iterator && it) = default;
template <typename P, typename = std::enable_if_t<not is_tensor<P>::value>>
iterator(P * data) : parent(data) {}
template <typename UP_P,
typename =
std::enable_if_t<is_tensor<typename UP_P::element_type>::value>>
iterator(UP_P && tensor) : parent(std::forward<UP_P>(tensor)) {}
iterator & operator=(const iterator & it) = default;
operator const_iterator<R>() {
return ConstConverterIteratorHelper<T, R>::convert(*this);
}
};
/* -------------------------------------------------------------------------- */
template <class T, bool is_scal>
template <typename R>
inline typename Array<T, is_scal>::template iterator<R>
Array<T, is_scal>::erase(const iterator<R> & it) {
T * curr = it.data();
UInt pos = (curr - values) / nb_component;
erase(pos);
iterator<R> rit = it;
return --rit;
}
#endif
} // namespace akantu
#endif /* __AKANTU_AKA_ARRAY_TMPL_HH__ */
diff --git a/src/mesh/mesh.hh b/src/mesh/mesh.hh
index fbe7d74ee..649b718d2 100644
--- a/src/mesh/mesh.hh
+++ b/src/mesh/mesh.hh
@@ -1,604 +1,603 @@
/**
* @file mesh.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Dana Christen <dana.christen@epfl.ch>
* @author David Simon Kammer <david.kammer@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Thu Jan 14 2016
*
* @brief the class representing the meshes
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MESH_HH__
#define __AKANTU_MESH_HH__
/* -------------------------------------------------------------------------- */
#include "aka_array.hh"
#include "aka_event_handler_manager.hh"
#include "aka_memory.hh"
#include "dumpable.hh"
#include "element.hh"
#include "element_class.hh"
#include "element_type_map.hh"
#include "group_manager.hh"
#include "mesh_data.hh"
#include "mesh_events.hh"
/* -------------------------------------------------------------------------- */
#include <set>
/* -------------------------------------------------------------------------- */
namespace akantu {
class Communicator;
class ElementSynchronizer;
class NodeSynchronizer;
} // namespace akantu
namespace akantu {
/* -------------------------------------------------------------------------- */
/* Mesh */
/* -------------------------------------------------------------------------- */
/**
* @class Mesh this contain the coordinates of the nodes in the Mesh.nodes
* Array, and the connectivity. The connectivity are stored in by element
* types.
*
* In order to loop on all element you have to loop on all types like this :
* @code{.cpp}
Mesh::type_iterator it = mesh.firstType(dim, ghost_type);
Mesh::type_iterator end = mesh.lastType(dim, ghost_type);
for(; it != end; ++it) {
UInt nb_element = mesh.getNbElement(*it);
const Array<UInt> & conn = mesh.getConnectivity(*it);
for(UInt e = 0; e < nb_element; ++e) {
...
}
}
@endcode
*/
class Mesh : protected Memory,
public EventHandlerManager<MeshEventHandler>,
public GroupManager,
public Dumpable {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
private:
/// default constructor used for chaining, the last parameter is just to
/// differentiate constructors
Mesh(UInt spatial_dimension, const ID & id, const MemoryID & memory_id,
Communicator & communicator);
public:
/// constructor that create nodes coordinates array
Mesh(UInt spatial_dimension, const ID & id = "mesh",
const MemoryID & memory_id = 0);
/// mesh not distributed and not using the default communicator
Mesh(UInt spatial_dimension, Communicator & communicator,
const ID & id = "mesh", const MemoryID & memory_id = 0);
/// constructor that use an existing nodes coordinates array, by knowing its
/// ID
// Mesh(UInt spatial_dimension, const ID & nodes_id, const ID & id,
// const MemoryID & memory_id = 0);
/**
* constructor that use an existing nodes coordinates
* array, by getting the vector of coordinates
*/
Mesh(UInt spatial_dimension, const std::shared_ptr<Array<Real>> & nodes,
const ID & id = "mesh", const MemoryID & memory_id = 0);
~Mesh() override;
/// @typedef ConnectivityTypeList list of the types present in a Mesh
using ConnectivityTypeList = std::set<ElementType>;
/// read the mesh from a file
void read(const std::string & filename,
const MeshIOType & mesh_io_type = _miot_auto);
/// write the mesh to a file
void write(const std::string & filename,
const MeshIOType & mesh_io_type = _miot_auto);
private:
/// initialize the connectivity to NULL and other stuff
void init();
/// function that computes the bounding box (fills xmin, xmax)
void computeBoundingBox();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// patitionate the mesh among the processors involved in their computation
virtual void distribute(Communicator & communicator);
virtual void distribute();
/// function to print the containt of the class
void printself(std::ostream & stream, int indent = 0) const override;
/// extract coordinates of nodes from an element
template <typename T>
inline void extractNodalValuesFromElement(const Array<T> & nodal_values,
T * elemental_values,
UInt * connectivity, UInt n_nodes,
UInt nb_degree_of_freedom) const;
// /// extract coordinates of nodes from a reversed element
// inline void extractNodalCoordinatesFromPBCElement(Real * local_coords,
// UInt * connectivity,
// UInt n_nodes);
- /// convert a element to a linearized element
- inline UInt elementToLinearized(const Element & elem) const;
-
- /// convert a linearized element to an element
- inline Element linearizedToElement(UInt linearized_element) const;
-
- /// update the types offsets array for the conversions
- // inline void updateTypesOffsets(const GhostType & ghost_type);
-
/// add a Array of connectivity for the type <type>.
inline void addConnectivityType(const ElementType & type,
const GhostType & ghost_type = _not_ghost);
/* ------------------------------------------------------------------------ */
template <class Event> inline void sendEvent(Event & event) {
// if(event.getList().size() != 0)
EventHandlerManager<MeshEventHandler>::sendEvent<Event>(event);
}
/// prepare the event to remove the elements listed
void eraseElements(const Array<Element> & elements);
/* ------------------------------------------------------------------------ */
template <typename T>
inline void removeNodesFromArray(Array<T> & vect,
const Array<UInt> & new_numbering);
/// initialize normals
void initNormals();
/// init facets' mesh
Mesh & initMeshFacets(const ID & id = "mesh_facets");
/// define parent mesh
void defineMeshParent(const Mesh & mesh);
/// get global connectivity array
void getGlobalConnectivity(ElementTypeMapArray<UInt> & global_connectivity);
public:
void getAssociatedElements(const Array<UInt> & node_list,
Array<Element> & elements);
private:
/// fills the nodes_to_elements structure
void fillNodesToElements();
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// get the id of the mesh
AKANTU_GET_MACRO(ID, Memory::id, const ID &);
/// get the id of the mesh
AKANTU_GET_MACRO(MemoryID, Memory::memory_id, const MemoryID &);
/// get the spatial dimension of the mesh = number of component of the
/// coordinates
AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt);
/// get the nodes Array aka coordinates
AKANTU_GET_MACRO(Nodes, *nodes, const Array<Real> &);
AKANTU_GET_MACRO_NOT_CONST(Nodes, *nodes, Array<Real> &);
/// get the normals for the elements
AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Normals, normals, Real);
/// get the number of nodes
AKANTU_GET_MACRO(NbNodes, nodes->size(), UInt);
/// get the Array of global ids of the nodes (only used in parallel)
AKANTU_GET_MACRO(GlobalNodesIds, *nodes_global_ids, const Array<UInt> &);
AKANTU_GET_MACRO_NOT_CONST(GlobalNodesIds, *nodes_global_ids, Array<UInt> &);
/// get the global id of a node
inline UInt getNodeGlobalId(UInt local_id) const;
/// get the global id of a node
inline UInt getNodeLocalId(UInt global_id) const;
/// get the global number of nodes
inline UInt getNbGlobalNodes() const;
/// get the nodes type Array
AKANTU_GET_MACRO(NodesType, nodes_type, const Array<NodeType> &);
protected:
AKANTU_GET_MACRO_NOT_CONST(NodesType, nodes_type, Array<NodeType> &);
public:
inline NodeType getNodeType(UInt local_id) const;
/// say if a node is a pure ghost node
inline bool isPureGhostNode(UInt n) const;
/// say if a node is pur local or master node
inline bool isLocalOrMasterNode(UInt n) const;
inline bool isLocalNode(UInt n) const;
inline bool isMasterNode(UInt n) const;
inline bool isSlaveNode(UInt n) const;
AKANTU_GET_MACRO(LowerBounds, lower_bounds, const Vector<Real> &);
AKANTU_GET_MACRO(UpperBounds, upper_bounds, const Vector<Real> &);
AKANTU_GET_MACRO(LocalLowerBounds, local_lower_bounds, const Vector<Real> &);
AKANTU_GET_MACRO(LocalUpperBounds, local_upper_bounds, const Vector<Real> &);
/// get the connectivity Array for a given type
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Connectivity, connectivities, UInt);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Connectivity, connectivities, UInt);
AKANTU_GET_MACRO(Connectivities, connectivities,
const ElementTypeMapArray<UInt> &);
/// get the number of element of a type in the mesh
inline UInt getNbElement(const ElementType & type,
const GhostType & ghost_type = _not_ghost) const;
/// get the number of element for a given ghost_type and a given dimension
inline UInt getNbElement(const UInt spatial_dimension = _all_dimensions,
const GhostType & ghost_type = _not_ghost,
const ElementKind & kind = _ek_not_defined) const;
// /// get the connectivity list either for the elements or the ghost elements
// inline const ConnectivityTypeList &
// getConnectivityTypeList(const GhostType & ghost_type = _not_ghost) const;
/// compute the barycenter of a given element
private:
inline void getBarycenter(UInt element, const ElementType & type,
Real * barycenter,
GhostType ghost_type = _not_ghost) const;
public:
inline void getBarycenter(const Element & element,
Vector<Real> & barycenter) const;
/// get the element connected to a subelement
const auto & getElementToSubelement() const;
/// get the element connected to a subelement
const auto &
getElementToSubelement(const ElementType & el_type,
const GhostType & ghost_type = _not_ghost) const;
/// get the element connected to a subelement
auto & getElementToSubelement(const ElementType & el_type,
const GhostType & ghost_type = _not_ghost);
/// get the subelement connected to an element
const auto &
getSubelementToElement(const ElementType & el_type,
const GhostType & ghost_type = _not_ghost) const;
/// get the subelement connected to an element
auto & getSubelementToElement(const ElementType & el_type,
const GhostType & ghost_type = _not_ghost);
/// register a new ElementalTypeMap in the MeshData
template <typename T>
inline ElementTypeMapArray<T> & registerData(const std::string & data_name);
template <typename T>
inline bool hasData(const ID & data_name, const ElementType & el_type,
const GhostType & ghost_type = _not_ghost) const;
/// get a name field associated to the mesh
template <typename T>
inline const Array<T> &
getData(const ID & data_name, const ElementType & el_type,
const GhostType & ghost_type = _not_ghost) const;
/// get a name field associated to the mesh
template <typename T>
inline Array<T> & getData(const ID & data_name, const ElementType & el_type,
const GhostType & ghost_type = _not_ghost);
/// tells if the data data_name is contained in the mesh or not
inline bool hasData(const ID & data_name) const;
/// get a name field associated to the mesh
template <typename T>
inline const ElementTypeMapArray<T> & getData(const ID & data_name) const;
/// get a name field associated to the mesh
template <typename T>
inline ElementTypeMapArray<T> & getData(const ID & data_name);
template <typename T>
ElementTypeMap<UInt> getNbDataPerElem(ElementTypeMapArray<T> & array,
const ElementKind & element_kind);
template <typename T>
dumper::Field * createFieldFromAttachedData(const std::string & field_id,
const std::string & group_name,
const ElementKind & element_kind);
/// templated getter returning the pointer to data in MeshData (modifiable)
template <typename T>
inline Array<T> &
getDataPointer(const std::string & data_name, const ElementType & el_type,
const GhostType & ghost_type = _not_ghost,
UInt nb_component = 1, bool size_to_nb_element = true,
bool resize_with_parent = false);
/// Facets mesh accessor
inline const Mesh & getMeshFacets() const;
inline Mesh & getMeshFacets();
/// Parent mesh accessor
inline const Mesh & getMeshParent() const;
inline bool isMeshFacets() const { return this->is_mesh_facets; }
/// defines is the mesh is distributed or not
inline bool isDistributed() const { return this->is_distributed; }
#ifndef SWIG
/// return the dumper from a group and and a dumper name
DumperIOHelper & getGroupDumper(const std::string & dumper_name,
const std::string & group_name);
#endif
/* ------------------------------------------------------------------------ */
/* Wrappers on ElementClass functions */
/* ------------------------------------------------------------------------ */
public:
/// get the number of nodes per element for a given element type
static inline UInt getNbNodesPerElement(const ElementType & type);
/// get the number of nodes per element for a given element type considered as
/// a first order element
static inline ElementType getP1ElementType(const ElementType & type);
/// get the kind of the element type
static inline ElementKind getKind(const ElementType & type);
/// get spatial dimension of a type of element
static inline UInt getSpatialDimension(const ElementType & type);
/// get number of facets of a given element type
static inline UInt getNbFacetsPerElement(const ElementType & type);
/// get number of facets of a given element type
static inline UInt getNbFacetsPerElement(const ElementType & type, UInt t);
+#ifndef SWIG
+
/// get local connectivity of a facet for a given facet type
static inline auto getFacetLocalConnectivity(const ElementType & type,
UInt t = 0);
/// get connectivity of facets for a given element
inline auto getFacetConnectivity(const Element & element, UInt t = 0) const;
+#endif
+
/// get the number of type of the surface element associated to a given
/// element type
static inline UInt getNbFacetTypes(const ElementType & type, UInt t = 0);
+#ifndef SWIG
+
/// get the type of the surface element associated to a given element
static inline constexpr auto getFacetType(const ElementType & type,
UInt t = 0);
/// get all the type of the surface element associated to a given element
static inline constexpr auto getAllFacetTypes(const ElementType & type);
+#endif
+
/// get the number of nodes in the given element list
static inline UInt getNbNodesPerElementList(const Array<Element> & elements);
/* ------------------------------------------------------------------------ */
/* Element type Iterator */
/* ------------------------------------------------------------------------ */
using type_iterator = ElementTypeMapArray<UInt, ElementType>::type_iterator;
using ElementTypesIteratorHelper =
ElementTypeMapArray<UInt, ElementType>::ElementTypesIteratorHelper;
template <typename... pack>
ElementTypesIteratorHelper elementTypes(pack &&... _pack) const;
inline type_iterator firstType(UInt dim = _all_dimensions,
GhostType ghost_type = _not_ghost,
ElementKind kind = _ek_regular) const {
return connectivities.firstType(dim, ghost_type, kind);
}
inline type_iterator lastType(UInt dim = _all_dimensions,
GhostType ghost_type = _not_ghost,
ElementKind kind = _ek_regular) const {
return connectivities.lastType(dim, ghost_type, kind);
}
AKANTU_GET_MACRO(ElementSynchronizer, *element_synchronizer,
const ElementSynchronizer &);
AKANTU_GET_MACRO_NOT_CONST(ElementSynchronizer, *element_synchronizer,
ElementSynchronizer &);
AKANTU_GET_MACRO(NodeSynchronizer, *node_synchronizer,
const NodeSynchronizer &);
AKANTU_GET_MACRO_NOT_CONST(NodeSynchronizer, *node_synchronizer,
NodeSynchronizer &);
// AKANTU_GET_MACRO_NOT_CONST(Communicator, *communicator, StaticCommunicator
// &);
#ifndef SWIG
AKANTU_GET_MACRO(Communicator, *communicator, const auto &);
AKANTU_GET_MACRO_NOT_CONST(Communicator, *communicator, auto &);
#endif
/* ------------------------------------------------------------------------ */
/* Private methods for friends */
/* ------------------------------------------------------------------------ */
private:
friend class MeshAccessor;
AKANTU_GET_MACRO(NodesPointer, *nodes, Array<Real> &);
/// get a pointer to the nodes_global_ids Array<UInt> and create it if
/// necessary
inline Array<UInt> & getNodesGlobalIdsPointer();
/// get a pointer to the nodes_type Array<Int> and create it if necessary
inline Array<NodeType> & getNodesTypePointer();
/// get a pointer to the connectivity Array for the given type and create it
/// if necessary
inline Array<UInt> &
getConnectivityPointer(const ElementType & type,
const GhostType & ghost_type = _not_ghost);
/// get the ghost element counter
inline Array<UInt> &
getGhostsCounters(const ElementType & type,
const GhostType & ghost_type = _ghost) {
AKANTU_DEBUG_ASSERT(ghost_type != _not_ghost,
"No ghost counter for _not_ghost elements");
return ghosts_counters(type, ghost_type);
}
/// get a pointer to the element_to_subelement Array for the given type and
/// create it if necessary
inline Array<std::vector<Element>> &
getElementToSubelementPointer(const ElementType & type,
const GhostType & ghost_type = _not_ghost);
/// get a pointer to the subelement_to_element Array for the given type and
/// create it if necessary
inline Array<Element> &
getSubelementToElementPointer(const ElementType & type,
const GhostType & ghost_type = _not_ghost);
AKANTU_GET_MACRO_NOT_CONST(MeshData, mesh_data, MeshData &);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
/// array of the nodes coordinates
std::shared_ptr<Array<Real>> nodes;
/// global node ids
std::unique_ptr<Array<UInt>> nodes_global_ids;
/// node type, -3 pure ghost, -2 master for the node, -1 normal node, i in
/// [0-N] slave node and master is proc i
Array<NodeType> nodes_type;
/// global number of nodes;
UInt nb_global_nodes{0};
/// all class of elements present in this mesh (for heterogenous meshes)
ElementTypeMapArray<UInt> connectivities;
/// count the references on ghost elements
ElementTypeMapArray<UInt> ghosts_counters;
/// map to normals for all class of elements present in this mesh
ElementTypeMapArray<Real> normals;
/// list of all existing types in the mesh
// ConnectivityTypeList type_set;
/// the spatial dimension of this mesh
UInt spatial_dimension{0};
/// types offsets
// Array<UInt> types_offsets;
// /// list of all existing types in the mesh
// ConnectivityTypeList ghost_type_set;
// /// ghost types offsets
// Array<UInt> ghost_types_offsets;
/// min of coordinates
Vector<Real> lower_bounds;
/// max of coordinates
Vector<Real> upper_bounds;
/// size covered by the mesh on each direction
Vector<Real> size;
/// local min of coordinates
Vector<Real> local_lower_bounds;
/// local max of coordinates
Vector<Real> local_upper_bounds;
/// Extra data loaded from the mesh file
MeshData mesh_data;
/// facets' mesh
std::unique_ptr<Mesh> mesh_facets;
/// parent mesh (this is set for mesh_facets meshes)
const Mesh * mesh_parent{nullptr};
/// defines if current mesh is mesh_facets or not
bool is_mesh_facets{false};
/// defines if the mesh is centralized or distributed
bool is_distributed{false};
/// Communicator on which mesh is distributed
Communicator * communicator;
/// Element synchronizer
std::unique_ptr<ElementSynchronizer> element_synchronizer;
/// Node synchronizer
std::unique_ptr<NodeSynchronizer> node_synchronizer;
using NodesToElements = std::vector<std::unique_ptr<std::set<Element>>>;
/// This info is stored to simplify the dynamic changes
NodesToElements nodes_to_elements;
};
/// standard output stream operator
inline std::ostream & operator<<(std::ostream & stream, const Mesh & _this) {
_this.printself(stream);
return stream;
}
} // namespace akantu
/* -------------------------------------------------------------------------- */
/* Inline functions */
/* -------------------------------------------------------------------------- */
#include "element_type_map_tmpl.hh"
#include "mesh_inline_impl.cc"
#endif /* __AKANTU_MESH_HH__ */
diff --git a/src/model/solid_mechanics/material.hh b/src/model/solid_mechanics/material.hh
index adb3385c7..19d09abde 100644
--- a/src/model/solid_mechanics/material.hh
+++ b/src/model/solid_mechanics/material.hh
@@ -1,682 +1,682 @@
/**
* @file material.hh
*
* @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Wed Nov 25 2015
*
* @brief Mother class for all materials
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_factory.hh"
#include "aka_voigthelper.hh"
#include "data_accessor.hh"
#include "integration_point.hh"
#include "parsable.hh"
#include "parser.hh"
/* -------------------------------------------------------------------------- */
#include "internal_field.hh"
#include "random_internal_field.hh"
/* -------------------------------------------------------------------------- */
#include "mesh_events.hh"
#include "solid_mechanics_model_event_handler.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MATERIAL_HH__
#define __AKANTU_MATERIAL_HH__
/* -------------------------------------------------------------------------- */
namespace akantu {
class Model;
class SolidMechanicsModel;
} // namespace akantu
namespace akantu {
/**
* Interface of all materials
* Prerequisites for a new material
* - inherit from this class
* - implement the following methods:
* \code
* virtual Real getStableTimeStep(Real h, const Element & element =
* ElementNull);
*
* virtual void computeStress(ElementType el_type,
* GhostType ghost_type = _not_ghost);
*
* virtual void computeTangentStiffness(const ElementType & el_type,
* Array<Real> & tangent_matrix,
* GhostType ghost_type = _not_ghost);
* \endcode
*
*/
class Material : public Memory,
- public DataAccessor<Element>,
+ public TemporarySwigDataAccessor<Element>,
public Parsable,
public MeshEventHandler,
protected SolidMechanicsModelEventHandler {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
#if __cplusplus > 199711L
Material(const Material & mat) = delete;
Material & operator=(const Material & mat) = delete;
#endif
/// Initialize material with defaults
Material(SolidMechanicsModel & model, const ID & id = "");
/// Initialize material with custom mesh & fe_engine
Material(SolidMechanicsModel & model, UInt dim, const Mesh & mesh,
FEEngine & fe_engine, const ID & id = "");
/// Destructor
~Material() override;
protected:
void initialize();
/* ------------------------------------------------------------------------ */
/* Function that materials can/should reimplement */
/* ------------------------------------------------------------------------ */
protected:
/// constitutive law
virtual void computeStress(__attribute__((unused)) ElementType el_type,
__attribute__((unused))
GhostType ghost_type = _not_ghost) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/// compute the tangent stiffness matrix
virtual void computeTangentModuli(__attribute__((unused))
const ElementType & el_type,
__attribute__((unused))
Array<Real> & tangent_matrix,
__attribute__((unused))
GhostType ghost_type = _not_ghost) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/// compute the potential energy
virtual void computePotentialEnergy(ElementType el_type,
GhostType ghost_type = _not_ghost);
/// compute the potential energy for an element
virtual void
computePotentialEnergyByElement(__attribute__((unused)) ElementType type,
__attribute__((unused)) UInt index,
__attribute__((unused))
Vector<Real> & epot_on_quad_points) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
virtual void updateEnergies(__attribute__((unused)) ElementType el_type,
__attribute__((unused))
GhostType ghost_type = _not_ghost) {}
virtual void updateEnergiesAfterDamage(__attribute__((unused))
ElementType el_type,
__attribute__((unused))
GhostType ghost_type = _not_ghost) {}
/// set the material to steady state (to be implemented for materials that
/// need it)
virtual void setToSteadyState(__attribute__((unused)) ElementType el_type,
__attribute__((unused))
GhostType ghost_type = _not_ghost) {}
/// function called to update the internal parameters when the modifiable
/// parameters are modified
virtual void updateInternalParameters() {}
public:
/// extrapolate internal values
virtual void extrapolateInternal(const ID & id, const Element & element,
const Matrix<Real> & points,
Matrix<Real> & extrapolated);
/// compute the p-wave speed in the material
virtual Real getPushWaveSpeed(__attribute__((unused))
const Element & element) const {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/// compute the s-wave speed in the material
virtual Real getShearWaveSpeed(__attribute__((unused))
const Element & element) const {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/// get a material celerity to compute the stable time step (default: is the
/// push wave speed)
virtual Real getCelerity(const Element & element) const {
return getPushWaveSpeed(element);
}
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
template <typename T>
void registerInternal(__attribute__((unused)) InternalField<T> & vect) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
template <typename T>
void unregisterInternal(__attribute__((unused)) InternalField<T> & vect) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/// initialize the material computed parameter
virtual void initMaterial();
/// compute the residual for this material
// virtual void updateResidual(GhostType ghost_type = _not_ghost);
/// assemble the residual for this material
virtual void assembleInternalForces(GhostType ghost_type);
/// save the stress in the previous_stress if needed
virtual void savePreviousState();
/// compute the stresses for this material
virtual void computeAllStresses(GhostType ghost_type = _not_ghost);
// virtual void
// computeAllStressesFromTangentModuli(GhostType ghost_type = _not_ghost);
virtual void computeAllCauchyStresses(GhostType ghost_type = _not_ghost);
/// set material to steady state
void setToSteadyState(GhostType ghost_type = _not_ghost);
/// compute the stiffness matrix
virtual void assembleStiffnessMatrix(GhostType ghost_type);
/// add an element to the local mesh filter
inline UInt addElement(const ElementType & type, UInt element,
const GhostType & ghost_type);
inline UInt addElement(const Element & element);
/// add many elements at once
void addElements(const Array<Element> & elements_to_add);
/// remove many element at once
void removeElements(const Array<Element> & elements_to_remove);
/// function to print the contain of the class
void printself(std::ostream & stream, int indent = 0) const override;
/**
* interpolate stress on given positions for each element by means
* of a geometrical interpolation on quadrature points
*/
void interpolateStress(ElementTypeMapArray<Real> & result,
const GhostType ghost_type = _not_ghost);
/**
* interpolate stress on given positions for each element by means
* of a geometrical interpolation on quadrature points and store the
* results per facet
*/
void interpolateStressOnFacets(ElementTypeMapArray<Real> & result,
ElementTypeMapArray<Real> & by_elem_result,
const GhostType ghost_type = _not_ghost);
/**
* function to initialize the elemental field interpolation
* function by inverting the quadrature points' coordinates
*/
void initElementalFieldInterpolation(
const ElementTypeMapArray<Real> & interpolation_points_coordinates);
/* ------------------------------------------------------------------------ */
/* Common part */
/* ------------------------------------------------------------------------ */
protected:
/* ------------------------------------------------------------------------ */
inline UInt getTangentStiffnessVoigtSize(UInt spatial_dimension) const;
/// compute the potential energy by element
void computePotentialEnergyByElements();
/// resize the intenals arrays
virtual void resizeInternals();
/* ------------------------------------------------------------------------ */
/* Finite deformation functions */
/* This functions area implementing what is described in the paper of Bathe */
/* et al, in IJNME, Finite Element Formulations For Large Deformation */
/* Dynamic Analysis, Vol 9, 353-386, 1975 */
/* ------------------------------------------------------------------------ */
protected:
/// assemble the residual
template <UInt dim> void assembleInternalForces(GhostType ghost_type);
/// Computation of Cauchy stress tensor in the case of finite deformation from
/// the 2nd Piola-Kirchhoff for a given element type
template <UInt dim>
void computeCauchyStress(ElementType el_type,
GhostType ghost_type = _not_ghost);
/// Computation the Cauchy stress the 2nd Piola-Kirchhoff and the deformation
/// gradient
template <UInt dim>
inline void computeCauchyStressOnQuad(const Matrix<Real> & F,
const Matrix<Real> & S,
Matrix<Real> & cauchy,
const Real & C33 = 1.0) const;
template <UInt dim>
void computeAllStressesFromTangentModuli(const ElementType & type,
GhostType ghost_type);
template <UInt dim>
void assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type);
/// assembling in finite deformation
template <UInt dim>
void assembleStiffnessMatrixNL(const ElementType & type,
GhostType ghost_type);
template <UInt dim>
void assembleStiffnessMatrixL2(const ElementType & type,
GhostType ghost_type);
/// Size of the Stress matrix for the case of finite deformation see: Bathe et
/// al, IJNME, Vol 9, 353-386, 1975
inline UInt getCauchyStressMatrixSize(UInt spatial_dimension) const;
/// Sets the stress matrix according to Bathe et al, IJNME, Vol 9, 353-386,
/// 1975
template <UInt dim>
inline void setCauchyStressMatrix(const Matrix<Real> & S_t,
Matrix<Real> & sigma);
/// write the stress tensor in the Voigt notation.
template <UInt dim>
inline void setCauchyStressArray(const Matrix<Real> & S_t,
Matrix<Real> & sigma_voight);
/* ------------------------------------------------------------------------ */
/* Conversion functions */
/* ------------------------------------------------------------------------ */
public:
template <UInt dim>
static inline void gradUToF(const Matrix<Real> & grad_u, Matrix<Real> & F);
static inline void rightCauchy(const Matrix<Real> & F, Matrix<Real> & C);
static inline void leftCauchy(const Matrix<Real> & F, Matrix<Real> & B);
template <UInt dim>
static inline void gradUToEpsilon(const Matrix<Real> & grad_u,
Matrix<Real> & epsilon);
template <UInt dim>
static inline void gradUToGreenStrain(const Matrix<Real> & grad_u,
Matrix<Real> & epsilon);
static inline Real stressToVonMises(const Matrix<Real> & stress);
protected:
/// converts global element to local element
inline Element convertToLocalElement(const Element & global_element) const;
/// converts local element to global element
inline Element convertToGlobalElement(const Element & local_element) const;
/// converts global quadrature point to local quadrature point
inline IntegrationPoint
convertToLocalPoint(const IntegrationPoint & global_point) const;
/// converts local quadrature point to global quadrature point
inline IntegrationPoint
convertToGlobalPoint(const IntegrationPoint & local_point) const;
/* ------------------------------------------------------------------------ */
/* DataAccessor inherited members */
/* ------------------------------------------------------------------------ */
public:
inline UInt getNbData(const Array<Element> & elements,
const SynchronizationTag & tag) const override;
inline void packData(CommunicationBuffer & buffer,
const Array<Element> & elements,
const SynchronizationTag & tag) const override;
inline void unpackData(CommunicationBuffer & buffer,
const Array<Element> & elements,
const SynchronizationTag & tag) override;
template <typename T>
inline void packElementDataHelper(const ElementTypeMapArray<T> & data_to_pack,
CommunicationBuffer & buffer,
const Array<Element> & elements,
const ID & fem_id = ID()) const;
template <typename T>
inline void unpackElementDataHelper(ElementTypeMapArray<T> & data_to_unpack,
CommunicationBuffer & buffer,
const Array<Element> & elements,
const ID & fem_id = ID());
/* ------------------------------------------------------------------------ */
/* MeshEventHandler inherited members */
/* ------------------------------------------------------------------------ */
public:
/* ------------------------------------------------------------------------ */
void onNodesAdded(const Array<UInt> &, const NewNodesEvent &) override{};
void onNodesRemoved(const Array<UInt> &, const Array<UInt> &,
const RemovedNodesEvent &) override{};
void onElementsAdded(const Array<Element> & element_list,
const NewElementsEvent & event) override;
void onElementsRemoved(const Array<Element> & element_list,
const ElementTypeMapArray<UInt> & new_numbering,
const RemovedElementsEvent & event) override;
void onElementsChanged(const Array<Element> &, const Array<Element> &,
const ElementTypeMapArray<UInt> &,
const ChangedElementsEvent &) override{};
/* ------------------------------------------------------------------------ */
/* SolidMechanicsModelEventHandler inherited members */
/* ------------------------------------------------------------------------ */
public:
virtual void beforeSolveStep();
virtual void afterSolveStep();
void onDamageIteration() override;
void onDamageUpdate() override;
void onDump() override;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
AKANTU_GET_MACRO(Name, name, const std::string &);
AKANTU_GET_MACRO(Model, model, const SolidMechanicsModel &)
AKANTU_GET_MACRO(ID, Memory::getID(), const ID &);
AKANTU_GET_MACRO(Rho, rho, Real);
AKANTU_SET_MACRO(Rho, rho, Real);
AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt);
/// return the potential energy for the subset of elements contained by the
/// material
Real getPotentialEnergy();
/// return the potential energy for the provided element
Real getPotentialEnergy(ElementType & type, UInt index);
/// return the energy (identified by id) for the subset of elements contained
/// by the material
virtual Real getEnergy(const std::string & energy_id);
/// return the energy (identified by id) for the provided element
virtual Real getEnergy(const std::string & energy_id, ElementType type,
UInt index);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementFilter, element_filter, UInt);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(GradU, gradu, Real);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(PotentialEnergy, potential_energy,
Real);
AKANTU_GET_MACRO(GradU, gradu, const ElementTypeMapArray<Real> &);
AKANTU_GET_MACRO(Stress, stress, const ElementTypeMapArray<Real> &);
AKANTU_GET_MACRO(ElementFilter, element_filter,
const ElementTypeMapArray<UInt> &);
AKANTU_GET_MACRO(FEEngine, fem, FEEngine &);
bool isNonLocal() const { return is_non_local; }
template <typename T>
const Array<T> & getArray(const ID & id, const ElementType & type,
const GhostType & ghost_type = _not_ghost) const;
template <typename T>
Array<T> & getArray(const ID & id, const ElementType & type,
const GhostType & ghost_type = _not_ghost);
template <typename T>
const InternalField<T> & getInternal(const ID & id) const;
template <typename T> InternalField<T> & getInternal(const ID & id);
template <typename T>
inline bool isInternal(const ID & id, const ElementKind & element_kind) const;
template <typename T>
ElementTypeMap<UInt>
getInternalDataPerElem(const ID & id, const ElementKind & element_kind) const;
bool isFiniteDeformation() const { return finite_deformation; }
bool isInelasticDeformation() const { return inelastic_deformation; }
template <typename T> inline void setParam(const ID & param, T value);
inline const Parameter & getParam(const ID & param) const;
template <typename T>
void flattenInternal(const std::string & field_id,
ElementTypeMapArray<T> & internal_flat,
const GhostType ghost_type = _not_ghost,
ElementKind element_kind = _ek_not_defined) const;
/// apply a constant eigengrad_u everywhere in the material
virtual void applyEigenGradU(const Matrix<Real> & prescribed_eigen_grad_u,
const GhostType = _not_ghost);
/// specify if the matrix need to be recomputed for this material
virtual bool hasStiffnessMatrixChanged() { return true; }
protected:
bool isInit() const { return is_init; }
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// boolean to know if the material has been initialized
bool is_init;
std::map<ID, InternalField<Real> *> internal_vectors_real;
std::map<ID, InternalField<UInt> *> internal_vectors_uint;
std::map<ID, InternalField<bool> *> internal_vectors_bool;
protected:
/// Link to the fem object in the model
FEEngine & fem;
/// Finite deformation
bool finite_deformation;
/// Finite deformation
bool inelastic_deformation;
/// material name
std::string name;
/// The model to witch the material belong
SolidMechanicsModel & model;
/// density : rho
Real rho;
/// spatial dimension
UInt spatial_dimension;
/// list of element handled by the material
ElementTypeMapArray<UInt> element_filter;
/// stresses arrays ordered by element types
InternalField<Real> stress;
/// eigengrad_u arrays ordered by element types
InternalField<Real> eigengradu;
/// grad_u arrays ordered by element types
InternalField<Real> gradu;
/// Green Lagrange strain (Finite deformation)
InternalField<Real> green_strain;
/// Second Piola-Kirchhoff stress tensor arrays ordered by element types
/// (Finite deformation)
InternalField<Real> piola_kirchhoff_2;
/// potential energy by element
InternalField<Real> potential_energy;
/// tell if using in non local mode or not
bool is_non_local;
/// tell if the material need the previous stress state
bool use_previous_stress;
/// tell if the material need the previous strain state
bool use_previous_gradu;
/// elemental field interpolation coordinates
InternalField<Real> interpolation_inverse_coordinates;
/// elemental field interpolation points
InternalField<Real> interpolation_points_matrices;
/// vector that contains the names of all the internals that need to
/// be transferred when material interfaces move
std::vector<ID> internals_to_transfer;
};
/// standard output stream operator
inline std::ostream & operator<<(std::ostream & stream,
const Material & _this) {
_this.printself(stream);
return stream;
}
} // namespace akantu
#include "material_inline_impl.cc"
#include "internal_field_tmpl.hh"
#include "random_internal_field_tmpl.hh"
/* -------------------------------------------------------------------------- */
/* Auto loop */
/* -------------------------------------------------------------------------- */
/// This can be used to automatically write the loop on quadrature points in
/// functions such as computeStress. This macro in addition to write the loop
/// provides two tensors (matrices) sigma and grad_u
#define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type) \
Array<Real>::matrix_iterator gradu_it = \
this->gradu(el_type, ghost_type) \
.begin(this->spatial_dimension, this->spatial_dimension); \
Array<Real>::matrix_iterator gradu_end = \
this->gradu(el_type, ghost_type) \
.end(this->spatial_dimension, this->spatial_dimension); \
\
this->stress(el_type, ghost_type) \
.resize(this->gradu(el_type, ghost_type).size()); \
\
Array<Real>::iterator<Matrix<Real>> stress_it = \
this->stress(el_type, ghost_type) \
.begin(this->spatial_dimension, this->spatial_dimension); \
\
if (this->isFiniteDeformation()) { \
this->piola_kirchhoff_2(el_type, ghost_type) \
.resize(this->gradu(el_type, ghost_type).size()); \
stress_it = this->piola_kirchhoff_2(el_type, ghost_type) \
.begin(this->spatial_dimension, this->spatial_dimension); \
} \
\
for (; gradu_it != gradu_end; ++gradu_it, ++stress_it) { \
Matrix<Real> & __attribute__((unused)) grad_u = *gradu_it; \
Matrix<Real> & __attribute__((unused)) sigma = *stress_it
#define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END }
/// This can be used to automatically write the loop on quadrature points in
/// functions such as computeTangentModuli. This macro in addition to write the
/// loop provides two tensors (matrices) sigma_tensor, grad_u, and a matrix
/// where the elemental tangent moduli should be stored in Voigt Notation
#define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_mat) \
Array<Real>::matrix_iterator gradu_it = \
this->gradu(el_type, ghost_type) \
.begin(this->spatial_dimension, this->spatial_dimension); \
Array<Real>::matrix_iterator gradu_end = \
this->gradu(el_type, ghost_type) \
.end(this->spatial_dimension, this->spatial_dimension); \
Array<Real>::matrix_iterator sigma_it = \
this->stress(el_type, ghost_type) \
.begin(this->spatial_dimension, this->spatial_dimension); \
\
tangent_mat.resize(this->gradu(el_type, ghost_type).size()); \
\
UInt tangent_size = \
this->getTangentStiffnessVoigtSize(this->spatial_dimension); \
Array<Real>::matrix_iterator tangent_it = \
tangent_mat.begin(tangent_size, tangent_size); \
\
for (; gradu_it != gradu_end; ++gradu_it, ++sigma_it, ++tangent_it) { \
Matrix<Real> & __attribute__((unused)) grad_u = *gradu_it; \
Matrix<Real> & __attribute__((unused)) sigma_tensor = *sigma_it; \
Matrix<Real> & tangent = *tangent_it
#define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END }
/* -------------------------------------------------------------------------- */
namespace akantu {
using MaterialFactory =
Factory<Material, ID, UInt, const ID &, SolidMechanicsModel &, const ID &>;
} // namespace akantu
#define INSTANTIATE_MATERIAL_ONLY(mat_name) \
template class mat_name<1>; \
template class mat_name<2>; \
template class mat_name<3>
#define MATERIAL_DEFAULT_PER_DIM_ALLOCATOR(id, mat_name) \
[](UInt dim, const ID &, SolidMechanicsModel & model, \
const ID & id) -> std::unique_ptr<Material> { \
switch (dim) { \
case 1: \
return std::make_unique<mat_name<1>>(model, id); \
case 2: \
return std::make_unique<mat_name<2>>(model, id); \
case 3: \
return std::make_unique<mat_name<3>>(model, id); \
default: \
AKANTU_EXCEPTION("The dimension " \
<< dim << "is not a valid dimension for the material " \
<< #id); \
} \
}
#define INSTANTIATE_MATERIAL(id, mat_name) \
INSTANTIATE_MATERIAL_ONLY(mat_name); \
static bool material_is_alocated_##id [[gnu::unused]] = \
MaterialFactory::getInstance().registerAllocator( \
#id, MATERIAL_DEFAULT_PER_DIM_ALLOCATOR(id, mat_name))
#endif /* __AKANTU_MATERIAL_HH__ */
diff --git a/src/model/solid_mechanics/materials/material_python/material_python.hh b/src/model/solid_mechanics/materials/material_python/material_python.hh
index badec2b76..9237401dd 100644
--- a/src/model/solid_mechanics/materials/material_python/material_python.hh
+++ b/src/model/solid_mechanics/materials/material_python/material_python.hh
@@ -1,120 +1,119 @@
/**
* @file material_python.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date creation: Fri Nov 13 2015
* @date last modification: Fri Nov 13 2015
*
* @brief Python material
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/**
* @file material_python.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
*/
/* -------------------------------------------------------------------------- */
#include "python_functor.hh"
/* -------------------------------------------------------------------------- */
#include "material.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MATERIAL_PYTHON_HH__
#define __AKANTU_MATERIAL_PYTHON_HH__
/* -------------------------------------------------------------------------- */
namespace akantu {
class MaterialPython : public Material, PythonFunctor {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
MaterialPython(SolidMechanicsModel & model, PyObject * obj,
const ID & id = "");
~MaterialPython() override = default;
- ;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
void registerInternals();
void initMaterial() override;
/// constitutive law for all element of a type
void computeStress(ElementType el_type,
GhostType ghost_type = _not_ghost) override;
/// constitutive law for a given quad point
template <typename it_type>
void computeStress(Matrix<Real> & grad_u, Matrix<Real> & sigma,
std::vector<it_type> & internal_iterators);
/// compute the tangent stiffness matrix for an element type
void computeTangentModuli(const ElementType & el_type,
Array<Real> & tangent_matrix,
GhostType ghost_type = _not_ghost) override;
/// compute the push wave speed of the material
Real getPushWaveSpeed(const Element & element) const override;
protected:
/// update the dissipated energy, must be called after the stress have been
/// computed
// virtual void updateEnergies(ElementType el_type, GhostType ghost_type){};
/// compute the tangent stiffness matrix for a given quadrature point
// inline void computeTangentModuliOnQuad(Matrix<Real> & tangent, Real &
// dam){};
/* ------------------------------------------------------------------------ */
/* DataAccessor inherited members */
/* ------------------------------------------------------------------------ */
public:
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
// virtual Real getEnergy(std::string type){};
// virtual Real getEnergy(std::string energy_id, ElementType type, UInt
// index){};
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
std::map<std::string, Real> local_params;
std::map<std::string, InternalField<Real> *> internals;
};
} // akantu
#endif /* __AKANTU_MATERIAL_PYTHON_HH__ */
diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh
index f015fb6d5..19bdd2dd2 100644
--- a/src/model/solid_mechanics/solid_mechanics_model.hh
+++ b/src/model/solid_mechanics/solid_mechanics_model.hh
@@ -1,562 +1,556 @@
/**
* @file solid_mechanics_model.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Tue Jul 27 2010
* @date last modification: Tue Jan 19 2016
*
* @brief Model of Solid Mechanics
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "boundary_condition.hh"
#include "data_accessor.hh"
#include "fe_engine.hh"
#include "model.hh"
#include "non_local_manager.hh"
#include "solid_mechanics_model_event_handler.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_SOLID_MECHANICS_MODEL_HH__
#define __AKANTU_SOLID_MECHANICS_MODEL_HH__
namespace akantu {
class Material;
class MaterialSelector;
class DumperIOHelper;
class NonLocalManager;
template <ElementKind kind, class IntegrationOrderFunctor>
class IntegratorGauss;
template <ElementKind kind> class ShapeLagrange;
} // namespace akantu
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
class SolidMechanicsModel
: public Model,
- public DataAccessor<Element>,
+ public TemporarySwigDataAccessor<Element>,
public DataAccessor<UInt>,
public BoundaryCondition<SolidMechanicsModel>,
public NonLocalManagerCallback,
public EventHandlerManager<SolidMechanicsModelEventHandler> {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
class NewMaterialElementsEvent : public NewElementsEvent {
public:
AKANTU_GET_MACRO_NOT_CONST(MaterialList, material, Array<UInt> &);
AKANTU_GET_MACRO(MaterialList, material, const Array<UInt> &);
protected:
Array<UInt> material;
};
using MyFEEngineType = FEEngineTemplate<IntegratorGauss, ShapeLagrange>;
protected:
using EventManager = EventHandlerManager<SolidMechanicsModelEventHandler>;
public:
SolidMechanicsModel(
Mesh & mesh, UInt spatial_dimension = _all_dimensions,
const ID & id = "solid_mechanics_model", const MemoryID & memory_id = 0,
const ModelType model_type = ModelType::_solid_mechanics_model);
~SolidMechanicsModel() override;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
protected:
/// initialize completely the model
void initFullImpl(
const ModelOptions & options = SolidMechanicsModelOptions()) override;
/// initialize all internal arrays for materials
virtual void initMaterials();
/// initialize the model
void initModel() override;
/// function to print the containt of the class
void printself(std::ostream & stream, int indent = 0) const override;
/// get some default values for derived classes
std::tuple<ID, TimeStepSolverType>
getDefaultSolverID(const AnalysisMethod & method) override;
/* ------------------------------------------------------------------------ */
/* Solver interface */
/* ------------------------------------------------------------------------ */
public:
/// assembles the stiffness matrix,
virtual void assembleStiffnessMatrix();
/// assembles the internal forces in the array internal_forces
virtual void assembleInternalForces();
protected:
/// callback for the solver, this adds f_{ext} - f_{int} to the residual
void assembleResidual() override;
/// get the type of matrix needed
MatrixType getMatrixType(const ID & matrix_id) override;
/// callback for the solver, this assembles different matrices
void assembleMatrix(const ID & matrix_id) override;
/// callback for the solver, this assembles the stiffness matrix
void assembleLumpedMatrix(const ID & matrix_id) override;
/// callback for the solver, this is called at beginning of solve
void predictor() override;
/// callback for the solver, this is called at end of solve
void corrector() override;
/// callback for the solver, this is called at beginning of solve
void beforeSolveStep() override;
/// callback for the solver, this is called at end of solve
void afterSolveStep() override;
/// Callback for the model to instantiate the matricees when needed
void initSolver(TimeStepSolverType time_step_solver_type,
NonLinearSolverType non_linear_solver_type) override;
protected:
/* ------------------------------------------------------------------------ */
TimeStepSolverType getDefaultSolverType() const override;
/* ------------------------------------------------------------------------ */
ModelSolverOptions
getDefaultSolverOptions(const TimeStepSolverType & type) const override;
public:
bool isDefaultSolverExplicit() {
return method == _explicit_lumped_mass ||
method == _explicit_consistent_mass;
}
protected:
/// update the current position vector
void updateCurrentPosition();
/* ------------------------------------------------------------------------ */
/* Materials (solid_mechanics_model_material.cc) */
/* ------------------------------------------------------------------------ */
public:
/// register an empty material of a given type
Material & registerNewMaterial(const ID & mat_name, const ID & mat_type,
const ID & opt_param);
/// reassigns materials depending on the material selector
virtual void reassignMaterial();
/// apply a constant eigen_grad_u on all quadrature points of a given material
virtual void applyEigenGradU(const Matrix<Real> & prescribed_eigen_grad_u,
const ID & material_name,
const GhostType ghost_type = _not_ghost);
protected:
/// register a material in the dynamic database
Material & registerNewMaterial(const ParserSection & mat_section);
/// read the material files to instantiate all the materials
void instantiateMaterials();
/// set the element_id_by_material and add the elements to the good materials
virtual void
assignMaterialToElements(const ElementTypeMapArray<UInt> * filter = nullptr);
/* ------------------------------------------------------------------------ */
/* Mass (solid_mechanics_model_mass.cc) */
/* ------------------------------------------------------------------------ */
public:
/// assemble the lumped mass matrix
void assembleMassLumped();
/// assemble the mass matrix for consistent mass resolutions
void assembleMass();
protected:
/// assemble the lumped mass matrix for local and ghost elements
void assembleMassLumped(GhostType ghost_type);
/// assemble the mass matrix for either _ghost or _not_ghost elements
void assembleMass(GhostType ghost_type);
/// fill a vector of rho
void computeRho(Array<Real> & rho, ElementType type, GhostType ghost_type);
/// compute the kinetic energy
Real getKineticEnergy();
Real getKineticEnergy(const ElementType & type, UInt index);
/// compute the external work (for impose displacement, the velocity should be
/// given too)
Real getExternalWork();
/* ------------------------------------------------------------------------ */
/* NonLocalManager inherited members */
/* ------------------------------------------------------------------------ */
protected:
void initializeNonLocal() override;
void updateDataForNonLocalCriterion(ElementTypeMapReal & criterion) override;
void computeNonLocalStresses(const GhostType & ghost_type) override;
void
insertIntegrationPointsInNeighborhoods(const GhostType & ghost_type) override;
/// update the values of the non local internal
void updateLocalInternal(ElementTypeMapReal & internal_flat,
const GhostType & ghost_type,
const ElementKind & kind) override;
/// copy the results of the averaging in the materials
void updateNonLocalInternal(ElementTypeMapReal & internal_flat,
const GhostType & ghost_type,
const ElementKind & kind) override;
/* ------------------------------------------------------------------------ */
/* Data Accessor inherited members */
/* ------------------------------------------------------------------------ */
public:
UInt getNbData(const Array<Element> & elements,
const SynchronizationTag & tag) const override;
void packData(CommunicationBuffer & buffer, const Array<Element> & elements,
const SynchronizationTag & tag) const override;
void unpackData(CommunicationBuffer & buffer, const Array<Element> & elements,
const SynchronizationTag & tag) override;
UInt getNbData(const Array<UInt> & dofs,
const SynchronizationTag & tag) const override;
void packData(CommunicationBuffer & buffer, const Array<UInt> & dofs,
const SynchronizationTag & tag) const override;
void unpackData(CommunicationBuffer & buffer, const Array<UInt> & dofs,
const SynchronizationTag & tag) override;
protected:
void
splitElementByMaterial(const Array<Element> & elements,
std::vector<Array<Element>> & elements_per_mat) const;
template <typename Operation>
void splitByMaterial(const Array<Element> & elements, Operation && op) const;
/* ------------------------------------------------------------------------ */
/* Mesh Event Handler inherited members */
/* ------------------------------------------------------------------------ */
protected:
void onNodesAdded(const Array<UInt> & nodes_list,
const NewNodesEvent & event) override;
void onNodesRemoved(const Array<UInt> & element_list,
const Array<UInt> & new_numbering,
const RemovedNodesEvent & event) override;
void onElementsAdded(const Array<Element> & nodes_list,
const NewElementsEvent & event) override;
void onElementsRemoved(const Array<Element> & element_list,
const ElementTypeMapArray<UInt> & new_numbering,
const RemovedElementsEvent & event) override;
void onElementsChanged(const Array<Element> &, const Array<Element> &,
const ElementTypeMapArray<UInt> &,
const ChangedElementsEvent &) override{};
/* ------------------------------------------------------------------------ */
/* Dumpable interface (kept for convenience) and dumper relative functions */
/* ------------------------------------------------------------------------ */
public:
virtual void onDump();
//! decide wether a field is a material internal or not
bool isInternal(const std::string & field_name,
const ElementKind & element_kind);
#ifndef SWIG
//! give the amount of data per element
virtual ElementTypeMap<UInt>
getInternalDataPerElem(const std::string & field_name,
const ElementKind & kind);
//! flatten a given material internal field
ElementTypeMapArray<Real> &
flattenInternal(const std::string & field_name, const ElementKind & kind,
const GhostType ghost_type = _not_ghost);
//! flatten all the registered material internals
void flattenAllRegisteredInternals(const ElementKind & kind);
#endif
dumper::Field * createNodalFieldReal(const std::string & field_name,
const std::string & group_name,
bool padding_flag) override;
dumper::Field * createNodalFieldBool(const std::string & field_name,
const std::string & group_name,
bool padding_flag) override;
dumper::Field * createElementalField(const std::string & field_name,
const std::string & group_name,
bool padding_flag,
const UInt & spatial_dimension,
const ElementKind & kind) override;
virtual void dump(const std::string & dumper_name);
virtual void dump(const std::string & dumper_name, UInt step);
virtual void dump(const std::string & dumper_name, Real time, UInt step);
void dump() override;
virtual void dump(UInt step);
virtual void dump(Real time, UInt step);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// return the dimension of the system space
AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt);
/// set the value of the time step
void setTimeStep(Real time_step, const ID & solver_id = "") override;
/// get the value of the conversion from forces/ mass to acceleration
AKANTU_GET_MACRO(F_M2A, f_m2a, Real);
/// set the value of the conversion from forces/ mass to acceleration
AKANTU_SET_MACRO(F_M2A, f_m2a, Real);
/// get the SolidMechanicsModel::displacement vector
AKANTU_GET_MACRO(Displacement, *displacement, Array<Real> &);
/// get the SolidMechanicsModel::previous_displacement vector
AKANTU_GET_MACRO(PreviousDisplacement, *previous_displacement, Array<Real> &);
/// get the SolidMechanicsModel::current_position vector \warn only consistent
/// after a call to SolidMechanicsModel::updateCurrentPosition
const Array<Real> & getCurrentPosition();
/// get the SolidMechanicsModel::increment vector \warn only consistent if
/// SolidMechanicsModel::setIncrementFlagOn has been called before
AKANTU_GET_MACRO(Increment, *displacement_increment, Array<Real> &);
/// get the lumped SolidMechanicsModel::mass vector
AKANTU_GET_MACRO(Mass, *mass, Array<Real> &);
/// get the SolidMechanicsModel::velocity vector
AKANTU_GET_MACRO(Velocity, *velocity, Array<Real> &);
/// get the SolidMechanicsModel::acceleration vector, updated by
/// SolidMechanicsModel::updateAcceleration
AKANTU_GET_MACRO(Acceleration, *acceleration, Array<Real> &);
/// get the SolidMechanicsModel::force vector (external forces)
AKANTU_GET_MACRO(Force, *external_force, Array<Real> &);
/// get the SolidMechanicsModel::internal_force vector (internal forces)
AKANTU_GET_MACRO(InternalForce, *internal_force, Array<Real> &);
/// get the SolidMechanicsModel::blocked_dofs vector
AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array<bool> &);
/// get the value of the SolidMechanicsModel::increment_flag
AKANTU_GET_MACRO(IncrementFlag, increment_flag, bool);
/// get a particular material (by material index)
inline Material & getMaterial(UInt mat_index);
/// get a particular material (by material index)
inline const Material & getMaterial(UInt mat_index) const;
/// get a particular material (by material name)
inline Material & getMaterial(const std::string & name);
/// get a particular material (by material name)
inline const Material & getMaterial(const std::string & name) const;
/// get a particular material id from is name
inline UInt getMaterialIndex(const std::string & name) const;
/// give the number of materials
inline UInt getNbMaterials() const { return materials.size(); }
/// give the material internal index from its id
Int getInternalIndexFromID(const ID & id) const;
/// compute the stable time step
Real getStableTimeStep();
/// get the energies
Real getEnergy(const std::string & energy_id);
/// compute the energy for energy
Real getEnergy(const std::string & energy_id, const ElementType & type,
UInt index);
AKANTU_GET_MACRO(MaterialByElement, material_index,
const ElementTypeMapArray<UInt> &);
/// vectors containing local material element index for each global element
/// index
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialByElement, material_index,
UInt);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialByElement, material_index, UInt);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialLocalNumbering,
material_local_numbering, UInt);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialLocalNumbering,
material_local_numbering, UInt);
AKANTU_GET_MACRO_NOT_CONST(MaterialSelector, *material_selector,
MaterialSelector &);
AKANTU_SET_MACRO(MaterialSelector, material_selector,
std::shared_ptr<MaterialSelector>);
/// Access the non_local_manager interface
AKANTU_GET_MACRO(NonLocalManager, *non_local_manager, NonLocalManager &);
/// get the FEEngine object to integrate or interpolate on the boundary
FEEngine & getFEEngineBoundary(const ID & name = "") override;
protected:
friend class Material;
protected:
/// compute the stable time step
Real getStableTimeStep(const GhostType & ghost_type);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// conversion coefficient form force/mass to acceleration
Real f_m2a;
/// displacements array
Array<Real> * displacement;
UInt displacement_release{0};
/// displacements array at the previous time step (used in finite deformation)
Array<Real> * previous_displacement{nullptr};
/// increment of displacement
Array<Real> * displacement_increment{nullptr};
/// lumped mass array
Array<Real> * mass{nullptr};
/// Check if materials need to recompute the mass array
bool need_to_reassemble_lumped_mass{true};
/// Check if materials need to recompute the mass matrix
bool need_to_reassemble_mass{true};
/// velocities array
Array<Real> * velocity{nullptr};
/// accelerations array
Array<Real> * acceleration{nullptr};
/// accelerations array
// Array<Real> * increment_acceleration;
/// external forces array
Array<Real> * external_force{nullptr};
/// internal forces array
Array<Real> * internal_force{nullptr};
/// array specifing if a degree of freedom is blocked or not
Array<bool> * blocked_dofs{nullptr};
/// array of current position used during update residual
Array<Real> * current_position{nullptr};
UInt current_position_release{0};
/// Arrays containing the material index for each element
ElementTypeMapArray<UInt> material_index;
/// Arrays containing the position in the element filter of the material
/// (material's local numbering)
ElementTypeMapArray<UInt> material_local_numbering;
-#ifdef SWIGPYTHON
-public:
-#endif
/// list of used materials
std::vector<std::unique_ptr<Material>> materials;
/// mapping between material name and material internal id
std::map<std::string, UInt> materials_names_to_id;
-#ifdef SWIGPYTHON
-protected:
-#endif
/// class defining of to choose a material
std::shared_ptr<MaterialSelector> material_selector;
/// flag defining if the increment must be computed or not
bool increment_flag;
/// tells if the material are instantiated
bool are_materials_instantiated;
using flatten_internal_map = std::map<std::pair<std::string, ElementKind>,
ElementTypeMapArray<Real> *>;
/// map a registered internals to be flattened for dump purposes
flatten_internal_map registered_internals;
/// non local manager
std::unique_ptr<NonLocalManager> non_local_manager;
};
/* -------------------------------------------------------------------------- */
namespace BC {
namespace Neumann {
using FromStress = FromHigherDim;
using FromTraction = FromSameDim;
} // namespace Neumann
} // namespace BC
} // namespace akantu
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#include "material.hh"
#include "parser.hh"
#include "solid_mechanics_model_inline_impl.cc"
#include "solid_mechanics_model_tmpl.hh"
/* -------------------------------------------------------------------------- */
#endif /* __AKANTU_SOLID_MECHANICS_MODEL_HH__ */
diff --git a/src/synchronizer/data_accessor.hh b/src/synchronizer/data_accessor.hh
index cebcb66b3..577e3e551 100644
--- a/src/synchronizer/data_accessor.hh
+++ b/src/synchronizer/data_accessor.hh
@@ -1,279 +1,283 @@
/**
* @file data_accessor.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Wed Sep 01 2010
* @date last modification: Tue Dec 08 2015
*
* @brief Interface of accessors for pack_unpack system
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "communication_buffer.hh"
#include "element.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_DATA_ACCESSOR_HH__
#define __AKANTU_DATA_ACCESSOR_HH__
namespace akantu {
class FEEngine;
} // akantu
namespace akantu {
class DataAccessorBase {
public:
DataAccessorBase() = default;
virtual ~DataAccessorBase() = default;
};
template <class T> class DataAccessor : public virtual DataAccessorBase {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
DataAccessor() = default;
~DataAccessor() override = default;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/**
* @brief get the number of data to exchange for a given array of T
* (elements or dofs) and a given akantu::SynchronizationTag
*/
virtual UInt getNbData(const Array<T> & elements,
const SynchronizationTag & tag) const = 0;
/**
* @brief pack the data for a given array of T (elements or dofs) and a given
* akantu::SynchronizationTag
*/
virtual void packData(CommunicationBuffer & buffer, const Array<T> & element,
const SynchronizationTag & tag) const = 0;
/**
* @brief unpack the data for a given array of T (elements or dofs) and a
* given akantu::SynchronizationTag
*/
virtual void unpackData(CommunicationBuffer & buffer,
const Array<T> & element,
const SynchronizationTag & tag) = 0;
};
/* -------------------------------------------------------------------------- */
/* Specialization */
/* -------------------------------------------------------------------------- */
template <> class DataAccessor<Element> : public virtual DataAccessorBase {
public:
DataAccessor() = default;
~DataAccessor() override = default;
virtual UInt getNbData(const Array<Element> & elements,
const SynchronizationTag & tag) const = 0;
virtual void packData(CommunicationBuffer & buffer,
const Array<Element> & element,
const SynchronizationTag & tag) const = 0;
virtual void unpackData(CommunicationBuffer & buffer,
const Array<Element> & element,
const SynchronizationTag & tag) = 0;
/* ------------------------------------------------------------------------ */
public:
template <typename T, bool pack_helper>
static void
packUnpackNodalDataHelper(Array<T> & data, CommunicationBuffer & buffer,
const Array<Element> & elements, const Mesh & mesh);
/* ------------------------------------------------------------------------ */
template <typename T, bool pack_helper>
static void packUnpackElementalDataHelper(
ElementTypeMapArray<T> & data_to_pack, CommunicationBuffer & buffer,
const Array<Element> & element, bool per_quadrature_point_data,
const FEEngine & fem);
/* ------------------------------------------------------------------------ */
template <typename T>
static void
packNodalDataHelper(const Array<T> & data, CommunicationBuffer & buffer,
const Array<Element> & elements, const Mesh & mesh) {
packUnpackNodalDataHelper<T, true>(const_cast<Array<T> &>(data), buffer,
elements, mesh);
}
template <typename T>
static inline void
unpackNodalDataHelper(Array<T> & data, CommunicationBuffer & buffer,
const Array<Element> & elements, const Mesh & mesh) {
packUnpackNodalDataHelper<T, false>(data, buffer, elements, mesh);
}
/* ------------------------------------------------------------------------ */
template <typename T>
static inline void
packElementalDataHelper(const ElementTypeMapArray<T> & data_to_pack,
CommunicationBuffer & buffer,
const Array<Element> & elements,
bool per_quadrature_point, const FEEngine & fem) {
packUnpackElementalDataHelper<T, true>(
const_cast<ElementTypeMapArray<T> &>(data_to_pack), buffer, elements,
per_quadrature_point, fem);
}
template <typename T>
static inline void
unpackElementalDataHelper(ElementTypeMapArray<T> & data_to_unpack,
CommunicationBuffer & buffer,
const Array<Element> & elements,
bool per_quadrature_point, const FEEngine & fem) {
packUnpackElementalDataHelper<T, false>(data_to_unpack, buffer, elements,
per_quadrature_point, fem);
}
};
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
template <> class DataAccessor<UInt> : public virtual DataAccessorBase {
public:
DataAccessor() = default;
~DataAccessor() override = default;
virtual UInt getNbData(const Array<UInt> & elements,
const SynchronizationTag & tag) const = 0;
virtual void packData(CommunicationBuffer & buffer,
const Array<UInt> & element,
const SynchronizationTag & tag) const = 0;
virtual void unpackData(CommunicationBuffer & buffer,
const Array<UInt> & element,
const SynchronizationTag & tag) = 0;
/* ------------------------------------------------------------------------ */
public:
template <typename T, bool pack_helper>
static void packUnpackDOFDataHelper(Array<T> & data,
CommunicationBuffer & buffer,
const Array<UInt> & dofs);
template <typename T>
static inline void packDOFDataHelper(const Array<T> & data_to_pack,
CommunicationBuffer & buffer,
const Array<UInt> & dofs) {
packUnpackDOFDataHelper<T, true>(const_cast<Array<T> &>(data_to_pack),
buffer, dofs);
}
template <typename T>
static inline void unpackDOFDataHelper(Array<T> & data_to_unpack,
CommunicationBuffer & buffer,
const Array<UInt> & dofs) {
packUnpackDOFDataHelper<T, false>(data_to_unpack, buffer, dofs);
}
};
/* -------------------------------------------------------------------------- */
template <typename T> class AddOperation {
public:
inline T operator()(T & a, T & b) { return a + b; };
};
template <typename T> class IdentityOperation {
public:
inline T & operator()(T &, T & b) { return b; };
};
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
template <template <class> class Op, class T>
class ReduceUIntDataAccessor : public virtual DataAccessor<UInt> {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
ReduceUIntDataAccessor(Array<T> & data, const SynchronizationTag & tag)
: data(data), tag(tag) {}
~ReduceUIntDataAccessor() override = default;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/* ------------------------------------------------------------------------ */
UInt getNbData(const Array<UInt> & elements,
const SynchronizationTag & tag) const override {
if (tag != this->tag)
return 0;
Vector<T> tmp(data.getNbComponent());
return elements.size() * CommunicationBuffer::sizeInBuffer(tmp);
}
/* ------------------------------------------------------------------------ */
void packData(CommunicationBuffer & buffer, const Array<UInt> & elements,
const SynchronizationTag & tag) const override {
if (tag != this->tag)
return;
auto data_it = data.begin(data.getNbComponent());
for (auto el : elements) {
buffer << Vector<T>(data_it[el]);
}
}
/* ------------------------------------------------------------------------ */
void unpackData(CommunicationBuffer & buffer, const Array<UInt> & elements,
const SynchronizationTag & tag) override {
if (tag != this->tag)
return;
auto data_it = data.begin(data.getNbComponent());
for (auto el : elements) {
Vector<T> unpacked(data.getNbComponent());
Vector<T> vect(data_it[el]);
buffer >> unpacked;
vect = oper(vect, unpacked);
}
}
protected:
/// data to (un)pack
Array<T> & data;
/// Tag to consider
SynchronizationTag tag;
/// reduction operator
Op<Vector<T>> oper;
};
/* -------------------------------------------------------------------------- */
template <class T>
using SimpleUIntDataAccessor = ReduceUIntDataAccessor<IdentityOperation, T>;
+template <class T>
+using TemporarySwigDataAccessor = DataAccessor<T>;
+
+
} // akantu
#endif /* __AKANTU_DATA_ACCESSOR_HH__ */