diff --git a/CMakeLists.txt b/CMakeLists.txt
index e9445554e..840669a96 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,281 +1,323 @@
#===============================================================================
# @file CMakeLists.txt
# @author Nicolas Richart <nicolas.richart@epfl.ch>
# @date Fri Jun 11 09:46:59 2010
#
# @section LICENSE
#
# <insert lisence here>
#
# @section DESCRIPTION
#
#===============================================================================
#===============================================================================
# _ ,-, _
# ,--, /: :\/': :`\/: :\
# |`; ' `,' `.; `: | _ _
# | | | ' | |. | | | |
# | : | F E | A S | T++ || __ _| | ____ _ _ __ | |_ _ _
# | :. | : | : | : | \ / _` | |/ / _` | '_ \| __| | | |
# \__/: :.. : :.. | :.. | ) | (_| | < (_| | | | | |_| |_| |
# `---',\___/,\___/ /' \__,_|_|\_\__,_|_| |_|\__|\__,_|
# `==._ .. . /'
# `-::-'
#===============================================================================
#===============================================================================
# CMake Project
#===============================================================================
cmake_minimum_required(VERSION 2.6)
project(AKANTU)
enable_language(CXX)
#===============================================================================
# Misc.
#===============================================================================
set(AKANTU_CMAKE_DIR "${AKANTU_SOURCE_DIR}/cmake")
set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake")
set(BUILD_SHARED_LIBS ON CACHE BOOL "Build shared libraries.")
#===============================================================================
# Version Number
#===============================================================================
# AKANTU version number. An even minor number corresponds to releases.
set(AKANTU_MAJOR_VERSION 0)
set(AKANTU_MINOR_VERSION 1)
set(AKANTU_BUILD_VERSION 0)
set(AKANTU_VERSION
"${AKANTU_MAJOR_VERSION}.${AKANTU_MINOR_VERSION}.${AKANTU_BUILD_VERSION}"
)
# Append the library version information to the library target properties
if(NOT AKANTU_NO_LIBRARY_VERSION)
set(AKANTU_LIBRARY_PROPERTIES ${AKANTU_LIBRARY_PROPERTIES}
VERSION "${AKANTU_VERSION}"
SOVERSION "${AKANTU_MAJOR_VERSION}.${AKANTU_MINOR_VERSION}"
)
endif(NOT AKANTU_NO_LIBRARY_VERSION)
#===============================================================================
# Options
#===============================================================================
option(AKANTU_DEBUG "Compiles akantu with the debug messages" ON)
option(AKANTU_DOCUMENTATION "Build source documentation using Doxygen." OFF)
# compilation switch
option(AKANTU_BUILD_CONTACT "Build the contact algorithm" OFF)
# features
option(AKANTU_USE_IOHELPER "Add IOHelper support in akantu" OFF)
option(AKANTU_USE_MPI "Add MPI support in akantu" OFF)
option(AKANTU_USE_SCOTCH "Add Scotch support in akantu" OFF)
+option(AKANTU_USE_MUMPS "Add Mumps support in akantu" OFF)
option(AKANTU_USE_BLAS "Use BLAS for arithmetic operations" OFF)
+option(AKANTU_USE_EPSN "Use EPSN streering environment" OFF)
#===============================================================================
# Options handling
#===============================================================================
# Debug
if(NOT AKANTU_DEBUG)
add_definitions(-DAKANTU_NDEBUG)
endif(NOT AKANTU_DEBUG)
# IOHelper
if(AKANTU_USE_IOHELPER)
find_package(IOHelper REQUIRED)
if(IOHELPER_FOUND)
add_definitions(-DAKANTU_USE_IOHELPER)
set(AKANTU_EXTERNAL_LIB_INCLUDE_PATH ${AKANTU_EXTERNAL_LIB_INCLUDE_PATH}
${IOHELPER_INCLUDE_PATH}
)
set(AKANTU_EXTERNAL_LIBRARIES ${AKANTU_EXTERNAL_LIBRARIES}
${IOHELPER_LIBRARIES}
)
endif(IOHELPER_FOUND)
endif(AKANTU_USE_IOHELPER)
# MPI
set(AKANTU_MPI_ON FALSE)
if(AKANTU_USE_MPI)
find_package(MPI REQUIRED)
if(MPI_FOUND)
add_definitions(-DAKANTU_USE_MPI)
set(AKANTU_EXTERNAL_LIB_INCLUDE_PATH ${AKANTU_EXTERNAL_LIB_INCLUDE_PATH}
${MPI_INCLUDE_PATH}
)
set(AKANTU_EXTERNAL_LIBRARIES ${AKANTU_EXTERNAL_LIBRARIES}
${MPI_LIBRARY}
${MPI_EXTRA_LIBRARY}
)
set(AKANTU_MPI_ON TRUE)
endif(MPI_FOUND)
endif(AKANTU_USE_MPI)
-
# Scotch
set(AKANTU_SCOTCH_ON FALSE)
if(AKANTU_USE_SCOTCH)
find_package(Scotch REQUIRED)
if(SCOTCH_FOUND)
add_definitions(-DAKANTU_USE_SCOTCH)
set(AKANTU_EXTERNAL_LIB_INCLUDE_PATH ${AKANTU_EXTERNAL_LIB_INCLUDE_PATH}
${SCOTCH_INCLUDE_PATH}
)
set(AKANTU_EXTERNAL_LIBRARIES ${AKANTU_EXTERNAL_LIBRARIES}
${SCOTCH_LIBRARIES}
)
set(AKANTU_SCOTCH_ON TRUE)
endif(SCOTCH_FOUND)
endif(AKANTU_USE_SCOTCH)
+# MUMPS
+set(AKANTU_MUMPS_ON FALSE)
+if(AKANTU_USE_MUMPS)
+ find_package(Mumps REQUIRED)
+ if(MUMPS_FOUND)
+ add_definitions(-DAKANTU_USE_MUMPS)
+ set(AKANTU_EXTERNAL_LIB_INCLUDE_PATH ${AKANTU_EXTERNAL_LIB_INCLUDE_PATH}
+ ${MUMPS_INCLUDE_PATH}
+ )
+ set(AKANTU_EXTERNAL_LIBRARIES ${AKANTU_EXTERNAL_LIBRARIES}
+ ${MUMPS_LIBRARIES}
+ )
+ set(AKANTU_MUMPS_ON TRUE)
+ endif(MUMPS_FOUND)
+endif(AKANTU_USE_MUMPS)
+
# BLAS
set(AKANTU_BLAS_ON FALSE)
if(AKANTU_USE_BLAS)
enable_language(Fortran)
find_package(BLAS)
if(BLAS_FOUND)
add_definitions(-DAKANTU_USE_CBLAS)
set(AKANTU_EXTERNAL_LIBRARIES ${AKANTU_EXTERNAL_LIBRARIES}
${BLAS_LIBRARIES}
)
set(AKANTU_BLAS_ON TRUE)
endif(BLAS_FOUND)
endif(AKANTU_USE_BLAS)
+# EPSN
+set(AKANTU_EPSN_ON FALSE)
+if(AKANTU_USE_EPSN)
+ find_package(EPSN)
+ if(EPSN_FOUND)
+ add_definitions(-DAKANTU_USE_EPSN)
+ set(AKANTU_EXTERNAL_LIB_INCLUDE_PATH ${AKANTU_EXTERNAL_LIB_INCLUDE_PATH}
+ ${EPSN_INCLUDE_DIR}
+ )
+ set(AKANTU_EXTERNAL_LIBRARIES ${AKANTU_EXTERNAL_LIBRARIES}
+ ${EPSN_LIBRARIES}
+ )
+ set(AKANTU_EPSN_ON TRUE)
+ endif(EPSN_FOUND)
+endif(AKANTU_USE_EPSN)
+
#===============================================================================
# Library
#===============================================================================
set(AKANTU_COMMON_SRC
common/aka_common.cc
common/aka_error.cc
common/aka_extern.cc
common/aka_static_memory.cc
common/aka_memory.cc
common/aka_vector.cc
common/aka_math.cc
fem/mesh.cc
fem/fem.cc
fem/element_class.cc
# model/integration_scheme/central_difference.cc
model/solid_mechanics/solid_mechanics_model.cc
model/solid_mechanics/solid_mechanics_model_mass.cc
model/solid_mechanics/solid_mechanics_model_boundary.cc
model/solid_mechanics/solid_mechanics_model_material.cc
model/solid_mechanics/material.cc
model/solid_mechanics/materials/material_elastic.cc
mesh_utils/mesh_io.cc
mesh_utils/mesh_io/mesh_io_msh.cc
mesh_utils/mesh_partition.cc
mesh_utils/mesh_utils.cc
-# solver/sparse_matrix.cc
+ solver/sparse_matrix.cc
+ solver/solver.cc
synchronizer/ghost_synchronizer.cc
synchronizer/synchronizer.cc
synchronizer/communicator.cc
synchronizer/static_communicator.cc
)
set(AKANTU_CONTACT_SRC
model/solid_mechanics/contact.cc
model/solid_mechanics/contact_search.cc
model/solid_mechanics/contact_neighbor_structure.cc
model/solid_mechanics/contact/contact_2d_explicit.cc
model/solid_mechanics/contact/contact_search_2d_explicit.cc
model/solid_mechanics/contact/regular_grid_neighbor_structure.cc
model/solid_mechanics/contact/contact_search_3d_explicit.cc
model/solid_mechanics/contact/contact_3d_explicit.cc
model/solid_mechanics/contact/grid_2d_neighbor_structure.cc
)
if(AKANTU_BUILD_CONTACT)
set(AKANTU_COMMON_SRC
${AKANTU_COMMON_SRC}
${AKANTU_CONTACT_SRC})
endif(AKANTU_BUILD_CONTACT)
if(AKANTU_USE_MPI AND MPI_FOUND)
set(AKANTU_COMMON_SRC
${AKANTU_COMMON_SRC}
synchronizer/static_communicator_mpi.cc
)
endif(AKANTU_USE_MPI AND MPI_FOUND)
if(AKANTU_SCOTCH_ON)
set(AKANTU_COMMON_SRC
${AKANTU_COMMON_SRC}
mesh_utils/mesh_partition/mesh_partition_scotch.cc
)
endif(AKANTU_SCOTCH_ON)
+if(AKANTU_MUMPS_ON)
+ set(AKANTU_COMMON_SRC
+ ${AKANTU_COMMON_SRC}
+ solver/solver_mumps.cc
+ )
+endif(AKANTU_MUMPS_ON)
+
+
set(AKANTU_INCLUDE_DIRS
common
fem/
mesh_utils/
mesh_utils/mesh_io/
mesh_utils/mesh_partition/
model/
model/integration_scheme
model/solid_mechanics
model/solid_mechanics/materials
model/solid_mechanics/contact
synchronizer/
solver/
)
include_directories(
${AKANTU_INCLUDE_DIRS}
${AKANTU_EXTERNAL_LIB_INCLUDE_PATH}
)
add_library(akantu ${AKANTU_COMMON_SRC})
set_target_properties(akantu PROPERTIES ${AKANTU_LIBRARY_PROPERTIES})
#===========================================================================
# Documentation
#===========================================================================
if(AKANTU_DOCUMENTATION)
find_package(Doxygen REQUIRED)
if(DOXYGEN_FOUND)
set(DOXYGEN_WARNINGS NO)
set(DOXYGEN_QUIET YES)
if(CMAKE_VERBOSE_MAKEFILE)
set(DOXYGEN_WARNINGS YES)
set(DOXYGEN_QUIET NO)
endif(CMAKE_VERBOSE_MAKEFILE)
add_subdirectory(doc/doxygen)
endif(DOXYGEN_FOUND)
endif(AKANTU_DOCUMENTATION)
#===============================================================================
# Tests
#===============================================================================
option(AKANTU_TESTS "Activate tests" OFF)
if(AKANTU_TESTS)
find_package(GMSH REQUIRED)
add_subdirectory(test)
endif(AKANTU_TESTS)
#===============================================================================
# Examples
#===============================================================================
option(AKANTU_EXAMPLES "Activate examples" OFF)
if(AKANTU_EXAMPLES)
find_package(GMSH REQUIRED)
add_subdirectory(examples)
endif(AKANTU_EXAMPLES)
diff --git a/cmake/FindScotch.cmake b/cmake/FindScotch.cmake
index 62e8488cb..b1b3a11a6 100644
--- a/cmake/FindScotch.cmake
+++ b/cmake/FindScotch.cmake
@@ -1,47 +1,47 @@
#===============================================================================
# @file FindScotch.cmake
# @author Nicolas Richart <nicolas.richart@epfl.ch>
# @date Tue Aug 25 16:53:57 2010
#
# @brief The find_package file for Scotch
#
# @section LICENSE
#
# <insert license here>
#
#===============================================================================
#===============================================================================
set(SCOTCH_LIBRARY "NOTFOUND" CACHE INTERNAL "Cleared" FORCE)
find_library(SCOTCH_LIBRARY scotch
PATHS ${SCOTCH_DIR}
PATH_SUFFIXES src/libscotch lib
)
find_library(SCOTCH_LIBRARY_ERR scotcherr
PATHS ${SCOTCH_DIR}
PATH_SUFFIXES src/libscotch lib
)
find_path(SCOTCH_INCLUDE_PATH scotch.h
PATHS ${SCOTCH_DIR}
PATH_SUFFIXES include scotch src/libscotch
)
#===============================================================================
mark_as_advanced(SCOTCH_LIBRARY)
mark_as_advanced(SCOTCH_LIBRARY_ERR)
mark_as_advanced(SCOTCH_INCLUDE_PATH)
set(SCOTCH_LIBRARIES_ALL ${SCOTCH_LIBRARY} ${SCOTCH_LIBRARY_ERR})
set(SCOTCH_LIBRARIES ${SCOTCH_LIBRARIES_ALL} CACHE INTERNAL "Libraries for scotch" FORCE)
#===============================================================================
include(FindPackageHandleStandardArgs)
find_package_handle_standard_args(SCOTCH DEFAULT_MSG
SCOTCH_LIBRARY SCOTCH_LIBRARY_ERR SCOTCH_INCLUDE_PATH)
#===============================================================================
if(NOT SCOTCH_FOUND)
- set(SCOTCH_DIR "" CACHE PATH "Location of IOHelper source directory.")
+ set(SCOTCH_DIR "" CACHE PATH "Location of Scotch library.")
endif(NOT SCOTCH_FOUND)
diff --git a/common/aka_common.hh b/common/aka_common.hh
index 2e8eae6b6..9529bc313 100644
--- a/common/aka_common.hh
+++ b/common/aka_common.hh
@@ -1,275 +1,278 @@
/**
* @file aka_common.hh
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Fri Jun 11 09:48:06 2010
*
* @namespace akantu
*
* @section LICENSE
*
* \<insert license here\>
*
* @section DESCRIPTION
*
* All common things to be included in the projects files
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_COMMON_HH__
#define __AKANTU_COMMON_HH__
/* -------------------------------------------------------------------------- */
#include <cmath>
#include <cstdlib>
#include <cstring>
/* -------------------------------------------------------------------------- */
#include <iostream>
#include <iomanip>
#include <string>
#include <exception>
#include <vector>
#include <map>
#include <set>
#include <list>
#include <limits>
#include <algorithm>
/* -------------------------------------------------------------------------- */
#define __BEGIN_AKANTU__ namespace akantu {
#define __END_AKANTU__ }
/* -------------------------------------------------------------------------- */
#include "aka_error.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
/* Common types */
/* -------------------------------------------------------------------------- */
typedef double Real;
typedef unsigned int UInt;
typedef int Int;
typedef std::string ID;
#ifdef AKANTU_NDEBUG
static const Real REAL_INIT_VALUE = 0;
#else
static const Real REAL_INIT_VALUE = std::numeric_limits<Real>::quiet_NaN();
#endif
/* -------------------------------------------------------------------------- */
/* Memory types */
/* -------------------------------------------------------------------------- */
typedef UInt MemoryID;
typedef ID VectorID;
/* -------------------------------------------------------------------------- */
/* Mesh/FEM/Model types */
/* -------------------------------------------------------------------------- */
typedef ID MeshID;
typedef ID FEMID;
typedef ID ModelID;
typedef ID MaterialID;
typedef ID SparseMatrixID;
+typedef ID SolverID;
typedef UInt Surface;
/// @enum ElementType type of element potentially contained in a Mesh
enum ElementType {
- _not_defined = 0,
+ _not_defined = 0,
_segment_2 = 1, /// first order segment
_segment_3 = 2, /// second order segment
- _triangle_3 = 3, /// first order triangle
- _triangle_6 = 4, /// second order triangle
- _tetrahedron_4 = 5, /// first order tetrahedron
- _tetrahedron_10 = 6, /// second order tetrahedron @remark not implemented yet
+ _triangle_3 = 3, /// first order triangle
+ _triangle_6 = 4, /// second order triangle
+ _tetrahedron_4 = 5, /// first order tetrahedron
+ _tetrahedron_10 = 6, /// second order tetrahedron @remark not implemented yet
+ _quadrangle_4, /// first order quadrangle
_max_element_type,
_point /// point only for some algorithm to be generic like mesh partitioning
};
/// @enum MaterialType different materials implemented
enum MaterialType {
_elastic = 0,
_max_material_type
};
/// @enum BoundaryFunctionType type of function passed for boundary conditions
enum BoundaryFunctionType {
_bft_forces,
_bft_stress
};
/// @enum SparseMatrixType type of sparse matrix used
enum SparseMatrixType {
_unsymmetric,
_symmetric
};
/* -------------------------------------------------------------------------- */
/* Contact */
/* -------------------------------------------------------------------------- */
typedef ID ContactID;
typedef ID ContactSearchID;
typedef ID ContactNeighborStructureID;
enum ContactType {
_ct_not_defined = 0,
_ct_2d_expli = 1,
_ct_3d_expli = 2
};
enum ContactSearchType {
_cst_not_defined = 0,
_cst_2d_expli = 1,
_cst_3d_expli = 2
};
enum ContactNeighborStructureType {
_cnst_not_defined = 0,
_cnst_regular_grid = 1,
_cnst_2d_grid = 2
};
/* -------------------------------------------------------------------------- */
/* Ghosts handling */
/* -------------------------------------------------------------------------- */
typedef ID SynchronizerID;
/// @enum CommunicatorType type of communication method to use
enum CommunicatorType {
_communicator_mpi,
_communicator_dummy
};
/// @enum GhostSynchronizationTag type of synchronizations
enum GhostSynchronizationTag {
/// SolidMechanicsModel tags
_gst_smm_mass, /// synchronization of the SolidMechanicsModel.mass
_gst_smm_residual, /// synchronization of the SolidMechanicsModel.current_position
_gst_smm_boundary, /// synchronization of the boundary, forces, velocities and displacement
/// Test tag
_gst_test
};
/// @enum GhostType type of ghost
enum GhostType {
_not_ghost,
_ghost,
_casper // not used but a real cute ghost
};
/// @enum SynchronizerOperation reduce operation that the synchronizer can perform
enum SynchronizerOperation {
_so_sum,
_so_min,
_so_max,
_so_null
};
/* -------------------------------------------------------------------------- */
/* Global defines */
/* -------------------------------------------------------------------------- */
#define AKANTU_MIN_ALLOCATION 2000
#define AKANTU_INDENT " "
/* -------------------------------------------------------------------------- */
#define AKANTU_SET_MACRO(name, variable, type) \
inline void set##name (type variable) { \
this->variable = variable; \
}
#define AKANTU_GET_MACRO(name, variable, type) \
inline type get##name () const { \
return variable; \
}
#define AKANTU_GET_MACRO_NOT_CONST(name, variable, type) \
inline type get##name () { \
return variable; \
}
#define AKANTU_GET_MACRO_BY_ELEMENT_TYPE(name, variable, type) \
inline type get##name (const ::akantu::ElementType & el_type) const { \
AKANTU_DEBUG_IN(); \
AKANTU_DEBUG_ASSERT(variable[el_type] != NULL, \
"No " << #variable << " of type " \
<< el_type << " in " << id); \
AKANTU_DEBUG_OUT(); \
return *variable[el_type]; \
}
/* -------------------------------------------------------------------------- */
//! standard output stream operator for ElementType
inline std::ostream & operator <<(std::ostream & stream, ElementType type)
{
switch(type)
{
- case _segment_2 : stream << "segment_2" ; break;
- case _segment_3 : stream << "segment_3" ; break;
- case _triangle_3 : stream << "triangle_3" ; break;
- case _triangle_6 : stream << "triangle_6" ; break;
- case _tetrahedron_4 : stream << "tetrahedron_4"; break;
- case _tetrahedron_10 : stream << "tetrahedron_10"; break;
- case _not_defined : stream << "undefined" ; break;
- case _max_element_type : stream << "ElementType(" << (int) type << ")"; break;
- case _point : stream << "point"; break;
+ case _segment_2 : stream << "segment_2" ; break;
+ case _segment_3 : stream << "segment_3" ; break;
+ case _triangle_3 : stream << "triangle_3" ; break;
+ case _triangle_6 : stream << "triangle_6" ; break;
+ case _tetrahedron_4 : stream << "tetrahedron_4" ; break;
+ case _tetrahedron_10 : stream << "tetrahedron_10"; break;
+ case _quadrangle_4 : stream << "quadrangle_4" ; break;
+ case _not_defined : stream << "undefined" ; break;
+ case _max_element_type : stream << "ElementType(" << (int) type << ")"; break;
+ case _point : stream << "point"; break;
}
return stream;
}
/* -------------------------------------------------------------------------- */
void initialize(int * argc, char *** argv);
/* -------------------------------------------------------------------------- */
void finalize ();
/* -------------------------------------------------------------------------- */
/*
* For intel compiler annoying remark
*/
#if defined(__INTEL_COMPILER)
/// remark #981: operands are evaluated in unspecified order
#pragma warning ( disable : 981 )
/// remark #383: value copied to temporary, reference to temporary used
#pragma warning ( disable : 383 )
#endif //defined(__INTEL_COMPILER)
/* -------------------------------------------------------------------------- */
/* string manipulation */
/* -------------------------------------------------------------------------- */
inline void to_lower(std::string & str) {
std::transform(str.begin(),
str.end(),
str.begin(),
(int(*)(int))std::tolower);
}
/* -------------------------------------------------------------------------- */
inline void trim(std::string & to_trim) {
size_t first = to_trim.find_first_not_of(" \t");
if (first != std::string::npos) {
size_t last = to_trim.find_last_not_of(" \t");
to_trim = to_trim.substr(first, last - first + 1);
} else to_trim = "";
}
__END_AKANTU__
#endif /* __AKANTU_COMMON_HH__ */
diff --git a/common/aka_error.cc b/common/aka_error.cc
index 6cf025553..933872364 100644
--- a/common/aka_error.cc
+++ b/common/aka_error.cc
@@ -1,107 +1,112 @@
/**
* @file aka_error.cc
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Sun Sep 5 21:03:37 2010
*
* @brief
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -------------------------------------------------------------------------- */
#include <csignal>
#include <cerrno>
#include <execinfo.h>
#include <cxxabi.h>
/* -------------------------------------------------------------------------- */
#include "aka_error.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
namespace debug {
/* ------------------------------------------------------------------------ */
void setDebugLevel(const DebugLevel & level) {
_debug_level = level;
}
+ /* ------------------------------------------------------------------------ */
+ const DebugLevel & getDebugLevel() {
+ return _debug_level;
+ }
+
/* ------------------------------------------------------------------------ */
void setParallelContext(int rank, int size) {
std::stringstream sstr;
sstr << "[" << std::setfill(' ') << std::right << std::setw(3) << (rank + 1) << "/" << size << "] ";
_parallel_context = sstr.str();
}
/* ------------------------------------------------------------------------ */
void initSignalHandler() {
struct sigaction action;
action.sa_handler = &printBacktrace;
sigemptyset(&(action.sa_mask));
action.sa_flags = SA_RESETHAND;
sigaction(SIGSEGV, &action, NULL);
}
/* ------------------------------------------------------------------------ */
static std::string demangle(const char* symbol) {
std::string trace(symbol);
std::string::size_type begin = trace.find_first_of('(') + 1;
std::string::size_type end = trace.find_last_of('+');
if (begin != std::string::npos && end != std::string::npos) {
std::string sub_trace = trace.substr(begin, end - begin);
int status;
size_t size;
std::string result;
char * demangled_name;
if ((demangled_name = abi::__cxa_demangle(sub_trace.c_str(), NULL, &size, &status)) != NULL) {
result = demangled_name;
free(demangled_name);
}
std::stringstream result_sstr;
result_sstr << trace.substr(0, begin) << result << trace.substr(end);
return result_sstr.str();
}
return trace;
}
/* ------------------------------------------------------------------------ */
void printBacktrace(int sig) {
AKANTU_DEBUG_INFO("Caught signal " << sig << "!");
void *array[10];
size_t size;
char **strings;
size_t i;
size = backtrace (array, 10);
strings = backtrace_symbols (array, size);
std::cerr << "BACKTRACE : " << size - 1 << " stack frames." <<std::endl;
/// -1 to remove the call to the printBacktrace function
for (i = 1; i < size; i++)
std::cerr << " " << demangle(strings[i]) << std::endl;;
free (strings);
std::cerr << "END BACKTRACE" << std::endl;
// int * segfault = NULL;
// *segfault = 0;
}
}
__END_AKANTU__
diff --git a/common/aka_error.hh b/common/aka_error.hh
index f5e10013e..19791c243 100644
--- a/common/aka_error.hh
+++ b/common/aka_error.hh
@@ -1,218 +1,220 @@
/**
* @file aka_error.hh
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Mon Jun 14 11:43:22 2010
*
* @brief error management and internal exceptions
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_ERROR_HH__
#define __AKANTU_ERROR_HH__
/* -------------------------------------------------------------------------- */
#include <ostream>
#include <sstream>
#ifdef AKANTU_USE_MPI
#include <mpi.h>
#endif
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
enum DebugLevel {
dbl0 = 0,
dblError = 0,
dblAssert = 0,
dbl1 = 1,
dblException = 1,
dblCritical = 1,
dbl2 = 2,
dblMajor = 2,
dbl3 = 3,
dblCall = 3,
dblSecondary = 3,
dblHead = 3,
dbl4 = 4,
dblWarning = 4,
dbl5 = 5,
dblInfo = 5,
dbl6 = 6,
dblIn = 6,
dblOut = 6,
dbl7 = 7,
dbl8 = 8,
dblTrace = 8,
dbl9 = 9,
dblAccessory = 9,
dbl10 = 10,
dbl100 = 100,
dblDump = 100
};
/* -------------------------------------------------------------------------- */
namespace debug {
extern std::ostream & _akantu_cout;
extern std::ostream & _akantu_debug_cout;
extern DebugLevel _debug_level;
extern std::string _parallel_context;
void setDebugLevel(const DebugLevel & level);
+ const DebugLevel & getDebugLevel();
void setParallelContext(int rank, int size);
void initSignalHandler();
void printBacktrace(int sig);
-}
-/* -------------------------------------------------------------------------- */
-/// exception class that can be thrown by akantu
-class Exception : public std::exception {
- /* ------------------------------------------------------------------------ */
- /* Constructors/Destructors */
- /* ------------------------------------------------------------------------ */
-public:
-
- //! full constructor
- Exception(std::string info, std::string file, unsigned int line) :
- _info(info), _file(file), _line(line) { }
-
- //! destructor
- virtual ~Exception() throw() {};
-
- /* ------------------------------------------------------------------------ */
- /* Methods */
- /* ------------------------------------------------------------------------ */
-public:
-
- virtual const char* what() const throw() {
- std::stringstream stream;
- stream << "akantu::Exception"
- << " : " << _info
- << " [" << _file << ":" << _line << "]";
- return stream.str().c_str();
- }
-
- virtual const char* info() const throw() {
- return _info.c_str();
- }
-
- /* ------------------------------------------------------------------------ */
- /* Class Members */
- /* ------------------------------------------------------------------------ */
-private:
-
- /// exception description and additionals
- std::string _info;
-
- /// file it is thrown from
- std::string _file;
-
- /// ligne it is thrown from
- unsigned int _line;
-};
+ /* -------------------------------------------------------------------------- */
+ /// exception class that can be thrown by akantu
+ class Exception : public std::exception {
+ /* ------------------------------------------------------------------------ */
+ /* Constructors/Destructors */
+ /* ------------------------------------------------------------------------ */
+ public:
+
+ //! full constructor
+ Exception(std::string info, std::string file, unsigned int line) :
+ _info(info), _file(file), _line(line) { }
+
+ //! destructor
+ virtual ~Exception() throw() {};
+
+ /* ------------------------------------------------------------------------ */
+ /* Methods */
+ /* ------------------------------------------------------------------------ */
+ public:
+
+ virtual const char* what() const throw() {
+ std::stringstream stream;
+ stream << "akantu::Exception"
+ << " : " << _info
+ << " [" << _file << ":" << _line << "]";
+ return stream.str().c_str();
+ }
+
+ virtual const char* info() const throw() {
+ return _info.c_str();
+ }
+
+ /* ------------------------------------------------------------------------ */
+ /* Class Members */
+ /* ------------------------------------------------------------------------ */
+ private:
+
+ /// exception description and additionals
+ std::string _info;
+
+ /// file it is thrown from
+ std::string _file;
+
+ /// ligne it is thrown from
+ unsigned int _line;
+ };
+};
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#define AKANTU_LOCATION "(" <<__FILE__ << ":" << __func__ << "():" << __LINE__ << ")"
#ifndef AKANTU_USE_MPI
#define AKANTU_EXIT(status) \
do { \
if (status != EXIT_SUCCESS) \
akantu::debug::printBacktrace(15); \
exit(status); \
} while(0)
#else
#define AKANTU_EXIT(status) \
do { \
if (status != EXIT_SUCCESS) \
akantu::debug::printBacktrace(15); \
MPI_Abort(MPI_COMM_WORLD, MPI_ERR_UNKNOWN); \
} while(0)
#endif
/* -------------------------------------------------------------------------- */
#define AKANTU_EXCEPTION(info) \
do { \
AKANTU_DEBUG(akantu::dblError, "!!! " << info); \
std::stringstream s_info; \
s_info << info ; \
- Exception ex(s_info.str(), __FILE__, __LINE__ ); \
+ ::akantu::debug::Exception ex(s_info.str(), __FILE__, __LINE__ ); \
throw ex; \
} while(0)
+
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_NDEBUG
#define AKANTU_DEBUG_TEST(level) (0)
#define AKANTU_DEBUG(level,info)
#define AKANTU_DEBUG_IN()
#define AKANTU_DEBUG_OUT()
#define AKANTU_DEBUG_INFO(info)
#define AKANTU_DEBUG_WARNING(info)
#define AKANTU_DEBUG_TRACE(info)
#define AKANTU_DEBUG_ASSERT(test,info)
#define AKANTU_DEBUG_ERROR(info) AKANTU_EXCEPTION(info)
/* -------------------------------------------------------------------------- */
#else
#define AKANTU_DEBUG(level,info) \
((::akantu::debug::_debug_level >= level) && \
(::akantu::debug::_akantu_debug_cout \
<< ::akantu::debug::_parallel_context \
<< info << " " \
<< AKANTU_LOCATION \
<< std::endl))
#define AKANTU_DEBUG_TEST(level) \
(::akantu::debug::_debug_level >= (level))
#define AKANTU_DEBUG_IN() \
AKANTU_DEBUG(::akantu::dblIn , "==> " << __func__ << "()")
#define AKANTU_DEBUG_OUT() \
AKANTU_DEBUG(::akantu::dblOut , "<== " << __func__ << "()")
#define AKANTU_DEBUG_INFO(info) \
AKANTU_DEBUG(::akantu::dblInfo , "--- " << info)
#define AKANTU_DEBUG_WARNING(info) \
AKANTU_DEBUG(::akantu::dblWarning, "??? " << info)
#define AKANTU_DEBUG_TRACE(info) \
AKANTU_DEBUG(::akantu::dblTrace , ">>> " << info)
#define AKANTU_DEBUG_ASSERT(test,info) \
do { \
if (!(test)) { \
AKANTU_DEBUG(::akantu::dblAssert, "assert [" << #test << "] " \
<< "!!! " << info); \
AKANTU_EXIT(EXIT_FAILURE); \
} \
} while(0)
#define AKANTU_DEBUG_ERROR(info) \
do { \
AKANTU_DEBUG(::akantu::dblError, "!!! " << info); \
AKANTU_EXIT(EXIT_FAILURE); \
} while(0)
#endif // AKANTU_NDEBUG
/* -------------------------------------------------------------------------- */
__END_AKANTU__
#endif /* __AKANTU_ERROR_HH__ */
// LocalWords: acessory
diff --git a/common/aka_extern.cc b/common/aka_extern.cc
index e95e27816..52c941ecd 100644
--- a/common/aka_extern.cc
+++ b/common/aka_extern.cc
@@ -1,48 +1,48 @@
/**
- * @file aka_extern.cpp
+ * @file aka_extern.cc
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Mon Jun 14 14:34:14 2010
*
* @brief initialisation of all global variables
* to insure the order of creation
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -------------------------------------------------------------------------- */
#include <ostream>
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
/** \todo write function to get this
* values from the environment or a config file
*/
/* -------------------------------------------------------------------------- */
/* error.hpp variables */
/* -------------------------------------------------------------------------- */
namespace debug {
/// standard output for debug messages
std::ostream & _akantu_debug_cout = std::cerr;
/// standard output for normal messages
std::ostream & _akantu_cout = std::cout;
/// debug level
DebugLevel _debug_level = (DebugLevel) 5;
/// parallel context used in debug messages
std::string _parallel_context = "";
}
/* -------------------------------------------------------------------------- */
__END_AKANTU__
diff --git a/common/aka_math.hh b/common/aka_math.hh
index 517d3268d..17a8ace94 100644
--- a/common/aka_math.hh
+++ b/common/aka_math.hh
@@ -1,163 +1,163 @@
/**
- * @file aka_math.h
+ * @file aka_math.hh
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Wed Jul 28 11:51:56 2010
*
* @brief mathematical operations
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_AKA_MATH_H__
#define __AKANTU_AKA_MATH_H__
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "aka_vector.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
class Math {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/* ------------------------------------------------------------------------ */
/* Matrix algebra */
/* ------------------------------------------------------------------------ */
/// @f$ y = A*x @f$
static void matrix_vector(UInt m, UInt n,
const Vector<Real> & A,
const Vector<Real> & x,
Vector<Real> & y);
/// @f$ y = A*x @f$
static inline void matrix_vector(UInt m, UInt n,
const Real * A,
const Real * x,
Real * y);
/// @f$ C = A*B @f$
static void matrix_matrix(UInt m, UInt n, UInt k,
const Vector<Real> & A,
const Vector<Real> & B,
Vector<Real> & C);
/// @f$ C = A*B @f$
static inline void matrix_matrix(UInt m, UInt n, UInt k,
const Real * A,
const Real * B,
Real * C);
/// @f$ C = A^t*B @f$
static inline void matrixt_matrix(UInt m, UInt n, UInt k,
const Real * A,
const Real * B,
Real * C);
/// @f$ C = A*B^t @f$
static inline void matrix_matrixt(UInt m, UInt n, UInt k,
const Real * A,
const Real * B,
Real * C);
/// @f$ C = A^t*B^t @f$
static inline void matrixt_matrixt(UInt m, UInt n, UInt k,
const Real * A,
const Real * B,
Real * C);
/// determinent of a 3x3 matrix
static inline Real det3(const Real * mat);
/// determinent of a 2x2 matrix
static inline Real det2(const Real * mat);
/// inverse a 3x3 matrix
static inline void inv3(const Real * mat, Real * inv);
/// inverse a 2x2 matrix
static inline void inv2(const Real * mat, Real * inv);
/// vector cross product
static inline void vectorProduct3(const Real * v1, const Real * v2, Real * res);
/// compute normal a normal to a vector
static inline void normal2(const Real * v1, Real * res);
/// compute normal a normal to a vector
static inline void normal3(const Real * v1,const Real * v2, Real * res);
/// normalize a vector
static inline void normalize2(Real * v);
/// normalize a vector
static inline void normalize3(Real * v);
/// return norm of a 2-vector
static inline Real norm2(const Real * v);
/// return norm of a 3-vector
static inline Real norm3(const Real * v);
/// return the dot product between 2 vectors in 2d
static inline Real vectorDot2(const Real * v1, const Real * v2);
/// return the dot product between 2 vectors in 3d
static inline Real vectorDot3(const Real * v1, const Real * v2);
/* ------------------------------------------------------------------------ */
/* Geometry */
/* ------------------------------------------------------------------------ */
/// distance in 2D between x and y
static inline Real distance_2d(const Real * x, const Real * y);
/// distance in 3D between x and y
static inline Real distance_3d(const Real * x, const Real * y);
/// radius of the in-circle of a triangle
static inline Real triangle_inradius(const Real * coord1, const Real * coord2, const Real * coord3);
/// radius of the in-circle of a tetrahedron
static inline Real tetrahedron_inradius(const Real * coord1, const Real * coord2, const Real * coord3, const Real * coord4);
/// volume of a tetrahedron
static inline Real tetrahedron_volume(const Real * coord1, const Real * coord2, const Real * coord3, const Real * coord4);
/// compute the barycenter of n points
static inline void barycenter(const Real * coord,
UInt nb_points, UInt spatial_dimension,
Real * barycenter);
/// vector between x and y
static inline void vector_2d(const Real * x, const Real * y, Real * vec);
/// vector pointing from x to y in 3 spatial dimension
static inline void vector_3d(const Real * x, const Real * y, Real * vec);
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#include "aka_math_inline_impl.cc"
__END_AKANTU__
#endif /* __AKANTU_AKA_MATH_H__ */
diff --git a/common/aka_memory.cc b/common/aka_memory.cc
index 3f6909c6d..ad54adc2c 100644
--- a/common/aka_memory.cc
+++ b/common/aka_memory.cc
@@ -1,46 +1,46 @@
/**
- * @file aka_memory.cpp
+ * @file aka_memory.cc
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Tue Jun 15 11:15:53 2010
*
* @brief static memory wrapper
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_memory.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
Memory::Memory(MemoryID memory_id) {
static_memory = StaticMemory::getStaticMemory();
this->memory_id = memory_id;
}
/* -------------------------------------------------------------------------- */
Memory::~Memory() {
if(StaticMemory::isInstantiated()) {
std::list<VectorID>::iterator it;
for (it = handeld_vectors_id.begin(); it != handeld_vectors_id.end(); ++it) {
AKANTU_DEBUG(dblAccessory, "Deleting the vector " << *it);
static_memory->sfree(memory_id, *it);
}
}
handeld_vectors_id.clear();
static_memory->destroy();
static_memory = NULL;
}
/* -------------------------------------------------------------------------- */
__END_AKANTU__
diff --git a/common/aka_static_memory.cc b/common/aka_static_memory.cc
index 12abe494b..9353d839b 100644
--- a/common/aka_static_memory.cc
+++ b/common/aka_static_memory.cc
@@ -1,138 +1,138 @@
/**
* @file aka_static_memory.cc
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Thu Jun 10 14:42:37 2010
*
* @brief Memory management
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -------------------------------------------------------------------------- */
#include <stdexcept>
#include <sstream>
/* -------------------------------------------------------------------------- */
#include "aka_static_memory.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
bool StaticMemory::is_instantiated = false;
StaticMemory * StaticMemory::single_static_memory = NULL;
UInt StaticMemory::nb_reference = 0;
/* -------------------------------------------------------------------------- */
StaticMemory * StaticMemory::getStaticMemory() {
if(!single_static_memory) {
single_static_memory = new StaticMemory();
is_instantiated = true;
}
nb_reference++;
return single_static_memory;
}
/* -------------------------------------------------------------------------- */
void StaticMemory::destroy() {
nb_reference--;
if(nb_reference == 0) {
delete single_static_memory;
}
}
/* -------------------------------------------------------------------------- */
StaticMemory::~StaticMemory() {
AKANTU_DEBUG_IN();
MemoryMap::iterator memory_it;
for(memory_it = memories.begin(); memory_it != memories.end(); ++memory_it) {
VectorMap::iterator vector_it;
for(vector_it = (memory_it->second).begin();
vector_it != (memory_it->second).end();
++vector_it) {
delete vector_it->second;
}
(memory_it->second).clear();
}
memories.clear();
is_instantiated = false;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void StaticMemory::sfree(const MemoryID & memory_id,
const VectorID & name) {
AKANTU_DEBUG_IN();
try {
VectorMap & vectors = const_cast<VectorMap &>(getMemory(memory_id));
VectorMap::iterator vector_it;
vector_it = vectors.find(name);
if(vector_it != vectors.end()) {
delete vector_it->second;
vectors.erase(vector_it);
AKANTU_DEBUG_INFO("Vector " << name << " removed from the static memory number " << memory_id);
AKANTU_DEBUG_OUT();
return;
}
- } catch (Exception e) {
+ } catch (debug::Exception e) {
AKANTU_DEBUG_INFO("The memory " << memory_id << " does not exist (perhaps already freed) ["
<< e.what() << "]");
AKANTU_DEBUG_OUT();
return;
}
AKANTU_DEBUG_INFO("The vector " << name << " does not exist (perhaps already freed)");
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void StaticMemory::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();
stream.precision(2);
stream << space << "StaticMemory [" << std::endl;
UInt nb_memories = memories.size();
stream << space << " + nb memories : " << nb_memories << std::endl;
Real tot_size = 0;
MemoryMap::const_iterator memory_it;
for(memory_it = memories.begin(); memory_it != memories.end(); ++memory_it) {
Real mem_size = 0;
stream << space << AKANTU_INDENT << "Memory [" << std::endl;
UInt mem_id = memory_it->first;
stream << space << AKANTU_INDENT << " + memory id : "
<< mem_id << std::endl;
UInt nb_vectors = (memory_it->second).size();
stream << space << AKANTU_INDENT << " + nb vectors : "
<< nb_vectors << std::endl;
stream.precision(prec);
VectorMap::const_iterator vector_it;
for(vector_it = (memory_it->second).begin();
vector_it != (memory_it->second).end();
++vector_it) {
(vector_it->second)->printself(stream, indent + 2);
mem_size += (vector_it->second)->getMemorySize() / 1024.;
}
stream.precision(2);
stream << space << AKANTU_INDENT << " + total size : " << mem_size << "kB" << std::endl;
stream << space << AKANTU_INDENT << "]" << std::endl;
tot_size += mem_size;
}
stream << space << " + total size : " << tot_size << "kB" << std::endl;
stream << space << "]" << std::endl;
stream.precision(prec);
}
/* -------------------------------------------------------------------------- */
__END_AKANTU__
diff --git a/common/aka_vector.cc b/common/aka_vector.cc
index 2d0a7430e..e93c86495 100644
--- a/common/aka_vector.cc
+++ b/common/aka_vector.cc
@@ -1,300 +1,309 @@
/**
* @file aka_vector.cc
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Thu Jun 17 15:14:24 2010
*
* @brief class vector
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "aka_vector.hh"
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
/* Functions VectorBase */
/* -------------------------------------------------------------------------- */
VectorBase::VectorBase(const VectorID & id) :
id(id), allocated_size(0), size(0), nb_component(1), size_of_type(0) {
}
/* -------------------------------------------------------------------------- */
VectorBase::~VectorBase() {
}
/* -------------------------------------------------------------------------- */
void VectorBase::printself(std::ostream & stream, int indent) const {
std::string space;
for(Int i = 0; i < indent; i++, space += AKANTU_INDENT);
stream << space << "VectorBase [" << std::endl;
stream << space << " + size : " << size << std::endl;
stream << space << " + nb component : " << nb_component << std::endl;
stream << space << " + allocated size : " << allocated_size << std::endl;
Real mem_size = (allocated_size * nb_component * size_of_type) / 1024.;
stream << space << " + size of type : " << size_of_type << "B" << std::endl;
stream << space << " + memory allocated : " << mem_size << "kB" << std::endl;
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
/* Functions Vector<T> */
/* -------------------------------------------------------------------------- */
template <class T> Vector<T>::Vector (UInt size,
UInt nb_component,
const VectorID & id) :
VectorBase(id), values(NULL) {
AKANTU_DEBUG_IN();
allocate(size, nb_component);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <class T> Vector<T>::Vector (UInt size,
UInt nb_component,
const T def_values[],
const VectorID & id) :
VectorBase(id), values(NULL) {
AKANTU_DEBUG_IN();
allocate(size, nb_component);
for (UInt i = 0; i < size; ++i) {
for (UInt j = 0; j < nb_component; ++j) {
values[i*nb_component + j] = def_values[j];
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <class T> Vector<T>::Vector (UInt size,
UInt nb_component,
const T & value,
const VectorID & id) :
VectorBase(id), values(NULL) {
AKANTU_DEBUG_IN();
allocate(size, nb_component);
for (UInt i = 0; i < nb_component*size; ++i) {
values[i] = value;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <class T> Vector<T>::Vector(const Vector<T>& vect, bool deep) {
AKANTU_DEBUG_IN();
this->id = vect.id;
if (deep) {
allocate(vect.size, vect.nb_component);
for (UInt i = 0; i < size; ++i) {
for (UInt j = 0; j < nb_component; ++j) {
values[i*nb_component + j] = vect.values[i*nb_component + j];
}
}
} else {
this->values = vect.values;
this->size = vect.size;
this->nb_component = vect.nb_component;
this->allocated_size = vect.allocated_size;
this->size_of_type = vect.size_of_type;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <class T> Vector<T>::~Vector () {
AKANTU_DEBUG_IN();
AKANTU_DEBUG(dblAccessory, "Freeing "
<< allocated_size*nb_component*sizeof(T) / 1024.
<< "kB (" << id <<")");
if(values)
free(values);
size = allocated_size = 0;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <class T> void Vector<T>::allocate(UInt size,
UInt nb_component) {
AKANTU_DEBUG_IN();
if (size == 0){
values = NULL;
} else {
values = static_cast<T*>(malloc(nb_component * size * sizeof(T)));
AKANTU_DEBUG_ASSERT(values != NULL,
"Cannot allocate "
<< nb_component * size * sizeof(T) / 1024.
<< "kB (" << id <<")");
}
if (values == NULL) {
this->size = this->allocated_size = 0;
} else {
AKANTU_DEBUG(dblAccessory, "Allocated "
<< size * nb_component * sizeof(T) / 1024.
<< "kB (" << id <<")");
this->size = this->allocated_size = size;
}
this->size_of_type = sizeof(T);
this->nb_component = nb_component;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <class T> void Vector<T>::resize(UInt new_size) {
AKANTU_DEBUG_IN();
/// free some memory
if(new_size <= allocated_size) {
if(allocated_size - new_size > AKANTU_MIN_ALLOCATION) {
AKANTU_DEBUG(dblAccessory, "Freeing "
<< (allocated_size - size)*nb_component*sizeof(T) / 1024.
<< "kB (" << id <<")");
/// Normally there are no allocation problem when reducing an array
T * tmp_ptr = static_cast<T*>(realloc(values, new_size * nb_component * sizeof(T)));
if(new_size != 0 && tmp_ptr == NULL) {
AKANTU_DEBUG_ERROR("Cannot free data (" << id << ")"
<< " [current allocated size : " << allocated_size << " | "
<< "requested size : " << new_size << "]");
}
values = tmp_ptr;
allocated_size = new_size;
}
size = new_size;
+
AKANTU_DEBUG_OUT();
return;
}
/// allocate more memory
UInt size_to_alloc = (new_size - allocated_size < AKANTU_MIN_ALLOCATION) ?
allocated_size + AKANTU_MIN_ALLOCATION : new_size;
T *tmp_ptr = static_cast<T*>(realloc(values, size_to_alloc * nb_component * sizeof(T)));
AKANTU_DEBUG_ASSERT(tmp_ptr != NULL,
"Cannot allocate "
<< size_to_alloc * nb_component * sizeof(T) / 1024.
<< "kB");
if (tmp_ptr == NULL) {
AKANTU_DEBUG_ERROR("Cannot allocate more data (" << id << ")"
<< " [current allocated size : " << allocated_size << " | "
<< "requested size : " << new_size << "]");
}
AKANTU_DEBUG(dblAccessory, "Allocating "
<< (size_to_alloc - allocated_size)*nb_component*sizeof(T) / 1024.
<< "kB");
allocated_size = size_to_alloc;
size = new_size;
values = tmp_ptr;
+
AKANTU_DEBUG_OUT();
}
+
/* -------------------------------------------------------------------------- */
-template <class T> void Vector<T>::extendComponentsInterlaced(UInt multiplicator,UInt stride) {
+template <class T> void Vector<T>::extendComponentsInterlaced(UInt multiplicator,
+ UInt stride) {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(multiplicator > 1,
"you can only extend a vector by changing the number of components");
AKANTU_DEBUG_ASSERT(nb_component%stride == 0,
"stride must divide actual number of components");
+
values = static_cast<T*>(realloc(values, nb_component*multiplicator*size* sizeof(T)));
+
UInt strided_component = nb_component/stride;
UInt new_component = strided_component * 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*strided_component/new_component;
for (UInt l = 0,p=stride-1; l < stride; ++l,--p) {
values[k*new_component*stride+m*stride+p] = values[k*strided_component*stride+n*stride+p];
}
}
}
+
nb_component = new_component*stride;
+
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <class T> Int Vector<T>::find(const T & elem) const {
AKANTU_DEBUG_IN();
UInt i = 0;
for (; (i < size) && (values[i] != elem); ++i);
AKANTU_DEBUG_OUT();
return (i == size) ? -1 : (Int) i;
}
/* -------------------------------------------------------------------------- */
template <> Int Vector<Real>::find(const Real & elem) const {
AKANTU_DEBUG_IN();
UInt i = 0;
Real epsilon = std::numeric_limits<Real>::epsilon();
for (; (i < size) && (fabs(values[i] - elem) <= epsilon); ++i);
AKANTU_DEBUG_OUT();
return (i == size) ? -1 : (Int) i;
}
/* -------------------------------------------------------------------------- */
template <class T> void Vector<T>::printself(std::ostream & stream, int indent) const {
std::string space;
for(Int i = 0; i < indent; i++, space += AKANTU_INDENT);
Real real_size = allocated_size * nb_component * size_of_type / 1024.0;
std::streamsize prec = stream.precision();
std::ios_base::fmtflags ff = stream.flags();
stream.setf (std::ios_base::showbase);
stream.precision(2);
stream << space << "Vector<" << typeid(T).name() << "> [" << std::endl;
stream << space << " + id : " << this->id << std::endl;
stream << space << " + size : " << this->size << std::endl;
stream << space << " + nb_component : " << this->nb_component << std::endl;
stream << space << " + allocated size : " << this->allocated_size << std::endl;
stream << space << " + memory size : "
<< real_size << "kB" << std::endl;
stream << space << " + address : " << std::hex << this->values
<< std::dec << std::endl;
stream.precision(prec);
stream.flags(ff);
if(AKANTU_DEBUG_TEST(dblDump)) {
stream << space << " + values : {";
for (UInt i = 0; i < this->size; ++i) {
stream << "{";
for (UInt j = 0; j < this->nb_component; ++j) {
stream << this->values[i*nb_component + j];
if(j != this->nb_component - 1) stream << ", ";
}
stream << "}";
if(i != this->size - 1) stream << ", ";
}
stream << "}" << std::endl;
}
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
class MaterialBase;
template class Vector<Int>;
template class Vector<UInt>;
template class Vector<Real>;
template class Vector<bool>;
__END_AKANTU__
diff --git a/common/aka_vector.hh b/common/aka_vector.hh
index 2c0986965..dece15897 100644
--- a/common/aka_vector.hh
+++ b/common/aka_vector.hh
@@ -1,195 +1,194 @@
/**
* @file aka_vector.hh
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Thu Jun 17 10:04:55 2010
*
* @brief class of vectors
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_VECTOR_HH__
#define __AKANTU_VECTOR_HH__
/* -------------------------------------------------------------------------- */
#include <typeinfo>
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
/// class that afford to store vectors in static memory
class VectorBase {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
VectorBase(const VectorID & id = "");
virtual ~VectorBase();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// get the amount of space allocated in bytes
inline UInt getMemorySize() const;
/// set the size to zero without freeing the allocated space
inline void empty();
/// function to print the containt of the class
virtual void printself(std::ostream & stream, int indent = 0) const;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
AKANTU_GET_MACRO(AllocatedSize, allocated_size, UInt);
AKANTU_GET_MACRO(Size, size, UInt);
AKANTU_GET_MACRO(NbComponent, nb_component, UInt);
AKANTU_GET_MACRO(ID, id, const VectorID &);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// id of the vector
VectorID id;
/// the size allocated
UInt allocated_size;
/// the size used
UInt size;
/// number of components
UInt nb_component;
/// size of the stored type
UInt size_of_type;
};
/* -------------------------------------------------------------------------- */
template<class T> class Vector : public VectorBase {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
/// Allocation of a new vector
Vector(UInt size = 0, UInt nb_component = 1,
const VectorID & id = "");
/// Allocation of a new vector with a default value
Vector(UInt size, UInt nb_component,
const T def_values[], const VectorID & id = "");
/// Allocation of a new vector with a default value
Vector(UInt size, UInt nb_component,
const T & value, const VectorID & id = "");
/// Copy constructor (deep copy if deep=true) \todo to implement
Vector(const Vector<T>& vect, bool deep = true);
virtual ~Vector();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// get jth componemt of the ith tuple in read-only
inline const T & get(UInt i, UInt j = 0) const;
/// get jth componemt of the ith tuple
inline T & at(UInt i, UInt j = 0);
/// add an element at the and of the vector with the value value for all
/// component
inline void push_back(const T& value);
/// add an element at the and of the vector
inline void push_back(const T new_elem[]);
/**
* remove an element and move the last one in the hole
* /!\ change the order in the vector
*/
inline void erase(UInt i);
/// change the size of the vector and allocate more memory if needed
void resize(UInt size);
/// change the number of components by interlacing data
void extendComponentsInterlaced(UInt multiplicator,UInt stride);
-
/// function to print the containt of the class
virtual void printself(std::ostream & stream, int indent = 0) const;
// Vector<T>& operator=(const Vector<T>& vect);
/// search elem in the vector, return the position of the first occurrence or
/// -1 if not found
Int find(const T & elem) const;
protected:
/// perform the allocation for the constructors
void allocate(UInt size, UInt nb_component = 1);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
UInt getSize() const{ return this->size; };
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
public:
/// array of values
T *values; // /!\ very dangerous
};
#include "aka_vector_inline_impl.cc"
/* -------------------------------------------------------------------------- */
/* Inline Functions Vector<T> */
/* -------------------------------------------------------------------------- */
template <class T>
inline std::ostream & operator<<(std::ostream & stream, const Vector<T> & _this)
{
_this.printself(stream);
return stream;
}
/* -------------------------------------------------------------------------- */
/* Inline Functions VectorBase */
/* -------------------------------------------------------------------------- */
inline std::ostream & operator<<(std::ostream & stream, const VectorBase & _this)
{
_this.printself(stream);
return stream;
}
__END_AKANTU__
#endif /* __AKANTU_VECTOR_HH__ */
diff --git a/doc/manual/figures/vectors.svg b/doc/manual/figures/vectors.svg
index 2c9ad1e05..2d413b4c8 100644
--- a/doc/manual/figures/vectors.svg
+++ b/doc/manual/figures/vectors.svg
@@ -1,264 +1,264 @@
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
<!-- Created with Inkscape (http://www.inkscape.org/) -->
<svg
xmlns:dc="http://purl.org/dc/elements/1.1/"
xmlns:cc="http://creativecommons.org/ns#"
xmlns:rdf="http://www.w3.org/1999/02/22-rdf-syntax-ns#"
xmlns:svg="http://www.w3.org/2000/svg"
xmlns="http://www.w3.org/2000/svg"
xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd"
xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape"
width="744.09448819"
height="1052.3622047"
id="svg2"
version="1.1"
inkscape:version="0.48.0 r9654"
- sodipodi:docname="New document 1">
+ sodipodi:docname="vectors.eps">
<defs
id="defs4">
<inkscape:path-effect
effect="skeletal"
id="path-effect3851"
is_visible="true"
pattern="m 6.4285714,23.571429 c 0,-2.76 2.24,-5 4.9999996,-5 2.76,0 5,2.24 5,5 0,2.76 -2.24,5 -5,5 -2.7599996,0 -4.9999996,-2.24 -4.9999996,-5 z"
copytype="single_stretched"
prop_scale="1"
scale_y_rel="false"
spacing="0"
normal_offset="0"
tang_offset="0"
prop_units="false"
vertical_pattern="false"
fuse_tolerance="0" />
</defs>
<sodipodi:namedview
id="base"
pagecolor="#ffffff"
bordercolor="#666666"
borderopacity="1.0"
inkscape:pageopacity="0.0"
inkscape:pageshadow="2"
inkscape:zoom="2.8"
inkscape:cx="203.77796"
inkscape:cy="969.78973"
inkscape:document-units="px"
inkscape:current-layer="layer1"
showgrid="true"
inkscape:snap-grids="true"
inkscape:snap-global="true"
inkscape:window-width="1680"
inkscape:window-height="973"
inkscape:window-x="1280"
inkscape:window-y="24"
inkscape:window-maximized="1">
<inkscape:grid
type="xygrid"
id="grid3883" />
</sodipodi:namedview>
<metadata
id="metadata7">
<rdf:RDF>
<cc:Work
rdf:about="">
<dc:format>image/svg+xml</dc:format>
<dc:type
rdf:resource="http://purl.org/dc/dcmitype/StillImage" />
<dc:title></dc:title>
</cc:Work>
</rdf:RDF>
</metadata>
<g
inkscape:label="Layer 1"
inkscape:groupmode="layer"
id="layer1">
<rect
style="fill:none;stroke:#000000;stroke-width:0.99999988;stroke-miterlimit:0;stroke-opacity:1;stroke-dasharray:none"
id="rect2985-0"
width="58.99984"
height="99"
x="50.500004"
y="52.862183"
ry="0" />
<path
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
d="m 50.500002,131.86218 58.999998,0"
id="path3027"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
d="m 50.500002,111.86218 58.999998,0"
id="path3027-0"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
d="m 50.500002,91.862183 58.999998,0"
id="path3027-9"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
d="m 50.500002,71.862183 58.999998,0"
id="path3027-6"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1;stroke-miterlimit:4;stroke-dasharray:2,2;stroke-dashoffset:0"
d="m 70.5,52.674248 0,99.187932"
id="path3061"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#000000;stroke-width:0.99999994;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1;stroke-miterlimit:4;stroke-dasharray:1.99999988,1.99999988;stroke-dashoffset:0"
d="m 90.5,52.674253 0,99.187927"
id="path3061-4"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#000000;stroke-width:0.30000001;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:none"
d="m 145.3,89.147897 0,-6.785714 59.7,0 0,7.5"
id="path3855"
inkscape:connector-curvature="0" />
<text
xml:space="preserve"
style="font-size:20px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="75"
y="40.362183"
id="text3857"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan3859"
x="75"
y="40.362183">c</tspan></text>
<text
xml:space="preserve"
style="font-size:20px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="23.714287"
y="107.6479"
id="text3861"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan3863"
x="23.714287"
y="107.6479">n</tspan></text>
<path
style="fill:none;stroke:#000000;stroke-width:0.30000001;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:none"
d="m 47.0091,152.19868 -6.7857,0.0135 -0.0734,-99.684358 7.49999,-0.01564"
id="path3855-9"
inkscape:connector-curvature="0" />
<rect
style="fill:none;stroke:#000000;stroke-width:1;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:none;stroke-dashoffset:0"
id="rect3887"
width="120"
height="20"
x="145"
y="92.362183" />
<path
style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:2, 2;stroke-dashoffset:0"
d="m 165,92.362183 0,19.999997"
id="path3891"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:2, 2;stroke-dashoffset:0"
d="m 185,92.362183 0,19.999997"
id="path3891-6"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#000000;stroke-width:1px;stroke-linecap:butt;stroke-linejoin:miter;stroke-opacity:1"
d="m 205,92.362183 0,19.999997"
id="path3891-3"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:2, 2;stroke-dashoffset:0"
d="m 225,92.362183 0,19.999997"
id="path3891-0"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:2, 2;stroke-dashoffset:0"
d="m 245,92.362183 0,19.999997"
id="path3891-1"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:2, 2;stroke-dashoffset:0"
d="m 320,92.362183 0,19.999997"
id="path3891-8"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#000000;stroke-width:1;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:2, 2;stroke-dashoffset:0"
d="m 340,92.362183 0,19.999997"
id="path3891-66"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#000000;stroke-width:0.30000001;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:none"
d="m 144.99999,80.933607 0,-13.571428 215.00001,0 0,15"
id="path3855-4"
inkscape:connector-curvature="0" />
<rect
style="fill:none;stroke:#000000;stroke-width:1;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:none;stroke-dashoffset:0"
id="rect3887-1"
width="60"
height="20"
x="300"
y="92.362183" />
<path
style="fill:none;stroke:#000000;stroke-width:3;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:none"
d="m 270,102.36218 5,0"
id="path3989"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#000000;stroke-width:3;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:none"
d="m 280,102.36218 5,0"
id="path3989-1"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#000000;stroke-width:3;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:none"
d="m 290,102.36218 5,0"
id="path3989-2"
inkscape:connector-curvature="0" />
<path
style="fill:none;stroke:#000000;stroke-width:0.30000001;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:4;stroke-opacity:1;stroke-dasharray:none"
d="m 50,49.147897 0,-6.785714 59.7,0 0,7.5"
id="path3855-4-0"
inkscape:connector-curvature="0" />
<text
xml:space="preserve"
style="font-size:20px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="170"
y="80.362183"
id="text3857-1"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan3859-2"
x="170"
y="80.362183">c</tspan></text>
<text
xml:space="preserve"
style="font-size:20px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="251"
y="65.362183"
id="text3861-0"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan3863-3"
x="251"
y="65.362183">n</tspan></text>
<text
xml:space="preserve"
style="font-size:14px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="70"
y="187.36218"
id="text4098"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan4100"
x="70"
y="187.36218">(a)</tspan></text>
<text
xml:space="preserve"
style="font-size:14px;font-style:normal;font-weight:normal;line-height:125%;letter-spacing:0px;word-spacing:0px;fill:#000000;fill-opacity:1;stroke:none;font-family:Sans"
x="235"
y="187.36218"
id="text4102"
sodipodi:linespacing="125%"><tspan
sodipodi:role="line"
id="tspan4104"
x="235"
y="187.36218">(b)</tspan></text>
</g>
</svg>
diff --git a/doc/manual/manual.tex b/doc/manual/manual.tex
index fb9627175..8f5d9106c 100644
--- a/doc/manual/manual.tex
+++ b/doc/manual/manual.tex
@@ -1,182 +1,202 @@
\documentclass[a4paper,11pt]{book}
\usepackage{hyperref}
\usepackage{xspace}
\usepackage[T1]{fontenc}
\usepackage{palatino}
% \usepackage[dvips,dvipdf]{graphics}
\usepackage{graphics}
\usepackage{epsfig}
\usepackage{url}
\usepackage{a4wide}
\usepackage{boxedminipage}
\usepackage{calc}
\usepackage{fancybox}
\usepackage{alltt}
\usepackage{moreverb}
\usepackage{longtable}
\usepackage{array}
\usepackage{listings}
\usepackage{color}
\usepackage{rotating}
\usepackage{sectionbox}
% \usepackage{ae}
% \usepackage{aecompl}
% \usepackage[cm]{aeguill}
% \usepackage{times}
% \usepackage{bookman}
\usepackage{palatino}
%\usepackage{verbatim}
\sloppy
\newcommand{\version}{0.1}
\newcommand{\code}[1]{{\tt{#1}}}
\title{\textbf{\Huge AKANTU}\\
\vspace{0.5cm}
\textbf{\huge User's Guide}\\
\vspace{1cm}
{\small \today{} --- Version \version}
}
\date{}
\begin{document}
\maketitle
\tableofcontents
\chapter{Data structures\label{chap:data-structure}}
\section{Vectors\label{sec:vectors}}
The Vector class is a template class that can store scalar types as Real, UInt,
-Int or bool. A Vector instance is defined by its size and its number of
-component. It also has an identifier and some extra internal variables for memory
-handling purpose.
+Int or bool. A Vector instance is defined by its size and its number of
+component. It also has an identifier and some extra internal variables for
+memory handling purpose.
\begin{itemize}
\item The size is the number of tupels stored in the Vector.
\item The number of component is the number of values stored for each tuple.
\end{itemize}
\begin{figure}[!htb]
\centering
\includegraphics{figures/vectors}
\caption{View of a Vector of size n and c compenents, (a) representation of
the vector, (b) representation of the storage of the same Vector.}
\label{fig:vectors}
\end{figure}
All the data are linerarized in memory in an array called values. This class
member is public. So for a Vector of size $n$ and $c$ component, to access the
$j^{th}$ component of th $i^{th}$ tuple, you have to get $values[i * c +
j]$. You can also access to this value with accessor \code{at(i, j)} or the
constant accessor \code{get(i, j)}.
If you want to store a matrix on each tuple you have to linearize it. For
exemple if you want to store a m * k matrix on each tuple you must specify c = m
* k. To access a particular value in a matrix of a tuple you will have to do
something like $values[ i * m * k + j_i * m + j_j]$.
\section{Vector storage convention within FE object\label{sec:FE-convention}}
-The point of this section is to describe the convention of storage for
-vectors intented to be passed to {\bf FE} object methods.
-Indeed a convention of necessary since gradient operators
-or integration loops will use vectors as input and ouput.
-Such vectors will be ordered with a specific convention that
-we intend to describe now.
-
-For vector objects, the size of the vector is always
-the number of nodes associated. The number of components
-is related to the order of the tensor considered. For a scalar
-field it is 1, for a vectorial field, the size of the vector
-is the number of components. For a $m\times n$ matrix field
-the number of components is $m\times n$.
-
-One common operation is the manipulation of continuum fields
-at the quadrature point positions. Here the size of the
-vector is $mn\_element \times nb\_quad\_points$ and the number
-of components is related to the stored object. For instance
-the method \verb$interpolateOnQuadraturePoints()$ take a nodal
-field stored in a vector($n\_nodes$,$dim$) and return a
-vector($n\_elem \times n\_quads$,$dim$).
-
-Basic gradient operations, like method \verb$gradientOnQuadraturePoints()$,
-will take as input a vector($n\_nodes$,$dim$) and return
-a vector($n\_elem \times n\_quads$,$dim\times spatial\_dimension$)
-where spatial dimension is the number of dimension in which the
-domain is embedded.
-
-In the same way the integration routines expect
-vector($n\_elem \times n\_quads$,$dim\times spatial\_dimension$)
-and will return vector($n\_elem$,$dim\times spatial\_dimension$).
-For non-Gaussian integrations, the input by be direction a nodal field.
-(At present time, only gaussian integrators are implemented within akantu).
-
-Last but not the least is the vectorial assembly process for which accept as input
-vector($n\_elem \times n\_quads \times n\_nodes\_per\_element$,$dim$)
+
+The point of this section is to describe the convention of storage for vectors
+intented to be passed to {\bf FE} object methods. Indeed a convention of
+necessary since gradient operators or integration loops will use vectors as
+input and ouput. Such vectors will be ordered with a specific convention that
+we intend to describe now.
+
+For vector objects, the size of the vector is always the number of nodes
+associated. The number of components is related to the order of the tensor
+considered. For a scalar field it is 1, for a vectorial field, the size of the
+vector is the number of components. For a $m\times n$ matrix field the number of
+components is $m\times n$.
+
+One common operation is the manipulation of continuum fields at the quadrature
+point positions. Here the size of the vector is $mn\_element \times
+nb\_quad\_points$ and the number of components is related to the stored
+object. For instance the method \verb$interpolateOnQuadraturePoints()$ take a
+nodal field stored in a vector($n\_nodes$,$dim$) and return a vector($n\_elem
+\times n\_quads$,$dim$).
+
+Basic gradient operations, like method \verb$gradientOnQuadraturePoints()$, will
+take as input a vector($n\_nodes$,$dim$) and return a vector($n\_elem \times
+n\_quads$,$dim \times spatial\_dimension$) where spatial dimension is the number
+of dimension in which the domain is embedded.
+
+In the same way the integration routines expect vector($n\_elem \times
+n\_quads$,$dim$) and will return vector($n\_elem$,$dim$). For non-Gaussian
+integrations, the input by be direction a nodal field. (At present time, only
+gaussian integrators are implemented within akantu).
+
+Last but not the least is the vectorial assembly process for which accept as
+input vector($n\_elem \times n\_quads \times n\_nodes\_per\_element$,$dim$)
and will return a vector($n\_nodes$,$dim$) nodal object.\\
-{\bf The general convention is that the number of component shall
-always be the size of the object manipulated at the lowest level.
-The object could be a per-element, of per quadrature point or even per node
-basis this will also apply as shown below. The figure \ref{vector-chain}
-shows the pattern of the vectors is the case a interpolation on quadrature points.}
+{\bf The general convention is that the number of component shall always be the
+ size of the object manipulated at the lowest level. The object could be a
+ per-element, of per quadrature point or even per node basis this will also
+ apply as shown below. The figures \ref{fig:vector-chain} and \ref{fig:interpolate-storage} shows the pattern of
+ the vectors is the case a interpolation on quadrature points.}
\begin{figure}
-\begin{equation}
-\left(
-\begin{array}{c}
-N_1 \\
-\vdots \\
-N_I
-\left\{
-\begin{array}{c}
-a_1 \\
-\vdots \\
-a_{dim} \\
-\end{array}
-\right. \\
-\vdots \\
-N_{nb\_nodes}
-\end{array}
-\right)
-\begin{array}{c}
-\Longrightarrow \\
-interpolate\\
-On\\
-Quadrature\\
-Points\\
-\end{array}
-\left(
-\begin{array}{c}
-E_1 \\
-\vdots \\
-E_e
-\left\{
-\begin{array}{c}
-q_1 \\
-\vdots \\
-\left.
-\begin{array}{c}
-a_1 \\
-\vdots \\
-a_{dim}
-\end{array}
-\right\} q_i \\
-\vdots \\
-q_p \\
-\end{array}
-\right. \\
-\vdots \\
-E_{nb\_elements}
-\end{array}
-\right)
-\end{equation}
-\caption{\label{vector-chain}Pattern of vectors manipulated during interpolation on quadrature points.
-Symbols N (resp. E, q) denotes nodes (resp. elements, quadrature points).}
-\end{figure}
+ \begin{equation}
+ \left(
+ \begin{array}{c}
+ N_1 \\
+ \vdots \\
+ N_I
+ \left\{
+ \begin{array}{c}
+ a_1 \\
+ \vdots \\
+ a_{dim} \\
+ \end{array}
+ \right. \\
+ \vdots \\
+ N_{nb\_nodes}
+ \end{array}
+ \right)
+ \begin{array}{c}
+ \Longrightarrow \\
+ interpolate\\
+ On\\
+ Quadrature\\
+ Points\\
+ \end{array}
+ \left(
+ \begin{array}{c}
+ E_1 \\
+ \vdots \\
+ E_e
+ \left\{
+ \begin{array}{c}
+ q_1 \\
+ \vdots \\
+ \left.
+ \begin{array}{c}
+ a_1 \\
+ \vdots \\
+ a_{dim}
+ \end{array}
+ \right\} q_i \\
+ \vdots \\
+ q_p \\
+ \end{array}
+ \right. \\
+ \vdots \\
+ E_{nb\_elements}
+ \end{array}
+ \right)
+ \end{equation}
+ \caption{\label{fig:vector-chain}Pattern of vectors manipulated during
+ interpolation on quadrature points. Symbols N (resp. E, q) denotes nodes
+ (resp. elements, quadrature points).}
+\end{figure}
+
+\begin{figure}[!htb]
+ \centering
+ \includegraphics{figures/interpolate}
+ \caption{Input and output vector of interpolateOnQuadraturePoints. The output
+ Vector has nb\_quadrature\_points tuples, the quadrature points are grouped by
+ elements (p is the number of quadrature points per element).}
+ \label{fig:interpolate-storage}
+\end{figure}
+
+\chapter{Models}
+
+\section{Solid Mechanics model}
+
+The solid mechanics model is a specifique implementation of the model
+interface. The model is instanciated for a given mesh. It will create is own FEM
+object to do the interpolation, gradient, integration and assemble
+operations. This model contain some Vectors that we will describe now.
+
+\begin{itemize}
+\item displacement
+\end{itemize}
+
\end{document}
diff --git a/fem/fem.hh b/fem/fem.hh
index 4035a2025..903a43a11 100644
--- a/fem/fem.hh
+++ b/fem/fem.hh
@@ -1,206 +1,208 @@
/**
* @file fem.hh
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Fri Jul 16 10:24:24 2010
*
* @brief FEM class
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_FEM_HH__
#define __AKANTU_FEM_HH__
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "aka_memory.hh"
#include "mesh.hh"
#include "element_class.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
class FEM : public Memory {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
FEM(Mesh & mesh, UInt spatial_dimension = 0,
FEMID id = "fem", MemoryID memory_id = 0);
virtual ~FEM();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// pre-compute all the shape functions, their derivatives and the jacobians
void initShapeFunctions(GhostType ghost_type = _not_ghost);
+ /// build the profile of the sparse matrix corresponding to the mesh
+ void initSparseMatrixProfile(SparseMatrixType sparse_matrix_type = _unsymmetric);
+
/// pre-compute normals on quadrature points
void computeNormalsOnQuadPoints(GhostType ghost_type = _not_ghost);
/// interpolate nodal values on the quadrature points
void interpolateOnQuadraturePoints(const Vector<Real> &u,
Vector<Real> &uq,
UInt nb_degre_of_freedom,
const ElementType & type,
GhostType ghost_type = _not_ghost,
const Vector<UInt> * filter_elements = NULL) const;
/// compute the gradient of u on the quadrature points
void gradientOnQuadraturePoints(const Vector<Real> &u,
Vector<Real> &nablauq,
UInt nb_degre_of_freedom,
const ElementType & type,
GhostType ghost_type = _not_ghost,
const Vector<UInt> * filter_elements = NULL) const;
/// integrate f for all elements of type "type"
void integrate(const Vector<Real> & f,
Vector<Real> &intf,
UInt nb_degre_of_freedom,
const ElementType & type,
GhostType ghost_type = _not_ghost,
const Vector<UInt> * filter_elements = NULL) const;
/// integrate f on the element "elem" of type "type"
inline void integrate(const Vector<Real> & f,
- Real *intf,
- UInt nb_degre_of_freedom,
- const ElementType & type,
- const UInt elem,
- GhostType ghost_type = _not_ghost) const;
+ Real *intf,
+ UInt nb_degre_of_freedom,
+ const Element & elem,
+ GhostType ghost_type = _not_ghost) const;
/// integrate a scalar value on all elements of type "type"
Real integrate(const Vector<Real> & f,
const ElementType & type,
GhostType ghost_type = _not_ghost,
const Vector<UInt> * filter_elements = NULL) const;
/// assemble vectors
void assembleVector(const Vector<Real> & elementary_vect,
Vector<Real> & nodal_values,
UInt nb_degre_of_freedom,
const ElementType & type,
GhostType ghost_type = _not_ghost,
const Vector<UInt> * filter_elements = NULL,
Real scale_factor = 1) const;
/// assemble matrix in the complete sparse matrix
void assembleMatrix() {};
/// function to print the containt of the class
virtual void printself(std::ostream & stream, int indent = 0) const;
private:
/// initialise the class
void init();
/// integrate on one element
inline void integrate(Real *f, Real *jac, Real * inte,
UInt nb_degre_of_freedom,
UInt nb_quadrature_points) const;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
AKANTU_GET_MACRO(ElementDimension, element_dimension, UInt);
/// get the mesh contained in the fem object
inline Mesh & getMesh() const;
/// get the size of the shapes returned by the element class
static inline UInt getShapeSize(const ElementType & type);
/// get the number of quadrature points
static inline UInt getNbQuadraturePoints(const ElementType & type);
/// get the size of the shapes derivatives returned by the element class
static inline UInt getShapeDerivativesSize(const ElementType & type);
/// get the in-radius of an element
static inline Real getElementInradius(Real * coord, const ElementType & type);
/// get a the shapes vector
AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Shapes, shapes, const Vector<Real> &);
/// get a the shapes derivatives vector
AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ShapesDerivatives, shapes_derivatives, const Vector<Real> &);
/// get the normals on quadrature points
AKANTU_GET_MACRO_BY_ELEMENT_TYPE(NormalsOnQuadPoints, normals_on_quad_points, const Vector<Real> &);
/// get a the ghost shapes vector
AKANTU_GET_MACRO_BY_ELEMENT_TYPE(GhostShapes, ghost_shapes,
const Vector<Real> &);
/// get a the ghost shapes derivatives vector
AKANTU_GET_MACRO_BY_ELEMENT_TYPE(GhostShapesDerivatives, ghost_shapes_derivatives,
const Vector<Real> &);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
/// id of the fem object
FEMID id;
/// spatial dimension of the problem
UInt element_dimension;
/// the mesh on which all computation are made
Mesh * mesh;
/// shape functions for all elements
ByElementTypeReal shapes;
/// shape derivatives for all elements
ByElementTypeReal shapes_derivatives;
/// jacobians for all elements
ByElementTypeReal jacobians;
/// normals at quadrature points
ByElementTypeReal normals_on_quad_points;
/// shape functions for all elements
ByElementTypeReal ghost_shapes;
/// shape derivatives for all elements
ByElementTypeReal ghost_shapes_derivatives;
/// jacobians for all elements
ByElementTypeReal ghost_jacobians;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#include "fem_inline_impl.cc"
/// standard output stream operator
inline std::ostream & operator <<(std::ostream & stream, const FEM & _this)
{
_this.printself(stream);
return stream;
}
__END_AKANTU__
#endif /* __AKANTU_FEM_HH__ */
diff --git a/mesh_utils/mesh_io.hh b/mesh_utils/mesh_io.hh
index 202e1e522..afd53342b 100644
--- a/mesh_utils/mesh_io.hh
+++ b/mesh_utils/mesh_io.hh
@@ -1,67 +1,69 @@
/**
* @file mesh_io.hh
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Fri Jun 18 10:27:42 2010
*
* @brief interface of a mesh io class, reader and writer
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MESH_IO_HH__
#define __AKANTU_MESH_IO_HH__
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "mesh.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
class MeshIO {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
MeshIO();
virtual ~MeshIO();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// read a mesh from the file
virtual void read(const std::string & filename, Mesh & mesh) = 0;
/// write a mesh to a file
virtual void write(const std::string & filename, const Mesh & mesh) = 0;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
bool canReadSurface;
bool canReadExtendedData;
// std::string filename;
// Mesh & mesh;
};
__END_AKANTU__
#endif /* __AKANTU_MESH_IO_HH__ */
+
+#include "mesh_io_msh.hh"
diff --git a/mesh_utils/mesh_io/mesh_io_msh.cc b/mesh_utils/mesh_io/mesh_io_msh.cc
index b7080961d..2e0667444 100644
--- a/mesh_utils/mesh_io/mesh_io_msh.cc
+++ b/mesh_utils/mesh_io/mesh_io_msh.cc
@@ -1,395 +1,396 @@
/**
* @file mesh_io_msh.cc
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Fri Jun 18 11:36:35 2010
*
* @brief Read/Write for MSH files generated by gmsh
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -----------------------------------------------------------------------------
Version (Legacy) 1.0
$NOD
number-of-nodes
node-number x-coord y-coord z-coord
...
$ENDNOD
$ELM
number-of-elements
elm-number elm-type reg-phys reg-elem number-of-nodes node-number-list
...
$ENDELM
-----------------------------------------------------------------------------
Version 2.1
$MeshFormat
version-number file-type data-size
$EndMeshFormat
$Nodes
number-of-nodes
node-number x-coord y-coord z-coord
...
$EndNodes
$Elements
number-of-elements
elm-number elm-type number-of-tags < tag > ... node-number-list
...
$EndElements
$PhysicalNames
number-of-names
physical-dimension physical-number "physical-name"
...
$EndPhysicalNames
$NodeData
number-of-string-tags
< "string-tag" >
...
number-of-real-tags
< real-tag >
...
number-of-integer-tags
< integer-tag >
...
node-number value ...
...
$EndNodeData
$ElementData
number-of-string-tags
< "string-tag" >
...
number-of-real-tags
< real-tag >
...
number-of-integer-tags
< integer-tag >
...
elm-number value ...
...
$EndElementData
$ElementNodeData
number-of-string-tags
< "string-tag" >
...
number-of-real-tags
< real-tag >
...
number-of-integer-tags
< integer-tag >
...
elm-number number-of-nodes-per-element value ...
...
$ElementEndNodeData
-----------------------------------------------------------------------------
elem-type
1: 2-node line.
2: 3-node triangle.
3: 4-node quadrangle.
4: 4-node tetrahedron.
5: 8-node hexahedron.
6: 6-node prism.
7: 5-node pyramid.
8: 3-node second order line
9: 6-node second order triangle
10: 9-node second order quadrangle
11: 10-node second order tetrahedron
12: 27-node second order hexahedron
13: 18-node second order prism
14: 14-node second order pyramid
15: 1-node point.
16: 8-node second order quadrangle
17: 20-node second order hexahedron
18: 15-node second order prism
19: 13-node second order pyramid
20: 9-node third order incomplete triangle
21: 10-node third order triangle
22: 12-node fourth order incomplete triangle
23: 15-node fourth order triangle
24: 15-node fifth order incomplete triangle
25: 21-node fifth order complete triangle
26: 4-node third order edge
27: 5-node fourth order edge
28: 6-node fifth order edge
29: 20-node third order tetrahedron
30: 35-node fourth order tetrahedron
31: 56-node fifth order tetrahedron
-------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#include <fstream>
/* -------------------------------------------------------------------------- */
#include "mesh_io_msh.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
/* Constant arrays initialisation */
/* -------------------------------------------------------------------------- */
UInt MeshIOMSH::_read_order[_max_element_type][MAX_NUMBER_OF_NODE_PER_ELEMENT] = {
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // _not_defined
{ 0, 1, 0, 0, 0, 0, 0, 0, 0, 0 }, // _line1
{ 0, 1, 2, 0, 0, 0, 0, 0, 0, 0 }, // _line2
{ 0, 1, 2, 0, 0, 0, 0, 0, 0, 0 }, // _triangle_3
{ 0, 1, 2, 3, 4, 5, 6, 0, 0, 0 }, // _triangle_6
{ 0, 1, 2, 3, 4, 0, 0, 0, 0, 0 }, // _tetrahedron_4
{ 0, 1, 2, 3, 4, 5, 6, 7, 9, 8 }, // _tetrahedron_10
+ { 0, 1, 2, 3, 0, 0, 0, 0, 0, 0 }, // _quadrangle_4
};
UInt MeshIOMSH::_msh_nodes_per_elem[16] =
{ 0, // element types began at 1
2, 3, 4, 4, 8, 6, 5, // 1st order
3, 6, 9, 10, 27, 18, 14, // 2st order
1
};
ElementType MeshIOMSH::_msh_to_akantu_element_types[16] =
{ _not_defined, // element types began at 1
- _segment_2, _triangle_3, _not_defined, _tetrahedron_4, // 1st order
+ _segment_2, _triangle_3, _quadrangle_4, _tetrahedron_4, // 1st order
_not_defined, _not_defined, _not_defined,
- _segment_3, _triangle_6, _not_defined, _tetrahedron_10, // 2nd order
+ _segment_3, _triangle_6, _not_defined, _tetrahedron_10, // 2nd order
_not_defined, _not_defined, _not_defined,
_not_defined
};
-MeshIOMSH::MSHElementType MeshIOMSH::_akantu_to_msh_element_types[_max_element_type
-] =
- { _msh_not_defined, // _not_defined
- _msh_segment_2, // _segment_2
- _msh_segment_3, // _segment_3
- _msh_triangle_3, // _triangle_3
- _msh_triangle_6, // _triangle_6
- _msh_tetrahedron_4, // _tetrahedron_4
- _msh_tetrahedron_10 // _tetrahedron_10
- };
+// MeshIOMSH::MSHElementType MeshIOMSH::_akantu_to_msh_element_types[_max_element_type] =
+// {
+// _msh_not_defined, // _not_defined
+// _msh_segment_2, // _segment_2
+// _msh_segment_3, // _segment_3
+// _msh_triangle_3, // _triangle_3
+// _msh_triangle_6, // _triangle_6
+// _msh_tetrahedron_4, // _tetrahedron_4
+// _msh_tetrahedron_10 // _tetrahedron_10
+// };
/* -------------------------------------------------------------------------- */
/* Methods Implentations */
/* -------------------------------------------------------------------------- */
MeshIOMSH::MeshIOMSH() {
canReadSurface = false;
canReadExtendedData = false;
}
/* -------------------------------------------------------------------------- */
MeshIOMSH::~MeshIOMSH() {
}
/* -------------------------------------------------------------------------- */
void MeshIOMSH::read(const std::string & filename, Mesh & mesh) {
std::ifstream infile;
infile.open(filename.c_str());
std::string line;
UInt first_node_number = 0, last_node_number = 0,
file_format = 1, current_line = 0;
if(!infile.good()) {
AKANTU_DEBUG_ERROR("Cannot open file " << filename);
}
while(infile.good()) {
std::getline(infile, line);
current_line++;
/// read the header
if(line == "$MeshFormat") {
std::getline(infile, line); /// the format line
std::getline(infile, line); /// the end of block line
current_line += 2;
file_format = 2;
}
/// read all nodes
if(line == "$Nodes" || line == "$NOD") {
UInt nb_nodes;
std::getline(infile, line);
std::stringstream sstr(line);
sstr >> nb_nodes;
current_line++;
Vector<Real> & nodes = const_cast<Vector<Real> &>(mesh.getNodes());
nodes.resize(nb_nodes);
UInt index;
Real coord[3];
UInt spatial_dimension = nodes.getNbComponent();
/// for each node, read the coordinates
for(UInt i = 0; i < nb_nodes; ++i) {
UInt offset = i * spatial_dimension;
std::getline(infile, line);
std::stringstream sstr_node(line);
sstr_node >> index >> coord[0] >> coord[1] >> coord[2];
current_line++;
first_node_number = first_node_number < index ? first_node_number : index;
last_node_number = last_node_number > index ? last_node_number : index;
/// read the coordinates
for(UInt j = 0; j < spatial_dimension; ++j)
nodes.values[offset + j] = coord[j];
}
std::getline(infile, line); /// the end of block line
}
/// read all elements
if(line == "$Elements" || line == "$ELM") {
UInt nb_elements;
std::getline(infile, line);
std::stringstream sstr(line);
sstr >> nb_elements;
current_line++;
Int index;
UInt msh_type;
ElementType akantu_type, akantu_type_old = _not_defined;
Vector<UInt> *connectivity = NULL;
UInt node_per_element = 0;
for(UInt i = 0; i < nb_elements; ++i) {
std::getline(infile, line);
std::stringstream sstr_elem(line);
current_line++;
sstr_elem >> index;
sstr_elem >> msh_type;
/// get the connectivity vector depending on the element type
akantu_type = _msh_to_akantu_element_types[msh_type];
if(akantu_type == _not_defined) continue;
if(akantu_type != akantu_type_old) {
connectivity = mesh.getConnectivityPointer(akantu_type);
connectivity->resize(0);
node_per_element = connectivity->getNbComponent();
akantu_type_old = akantu_type;
}
/// read tags informations
if(file_format == 2) {
UInt nb_tags;
sstr_elem >> nb_tags;
for(UInt j = 0; j < nb_tags; ++j) {
Int tag;
sstr_elem >> tag; ///@todo read to get extended information on elements
}
} else if (file_format == 1) {
Int tag;
sstr_elem >> tag; //reg-phys
sstr_elem >> tag; //reg-elem
sstr_elem >> tag; //number-of-nodes
}
UInt local_connect[node_per_element];
for(UInt j = 0; j < node_per_element; ++j) {
UInt node_index;
sstr_elem >> node_index;
AKANTU_DEBUG_ASSERT(node_index <= last_node_number,
"Node number not in range : line " << current_line);
node_index -= first_node_number + 1;
local_connect[_read_order[akantu_type][j]] = node_index;
}
connectivity->push_back(local_connect);
}
std::getline(infile, line); /// the end of block line
}
if((line[0] == '$') && (line.find("End") == std::string::npos)) {
AKANTU_DEBUG_WARNING("Unsuported block_kind " << line
<< " at line " << current_line);
}
}
infile.close();
}
/* -------------------------------------------------------------------------- */
void MeshIOMSH::write(const std::string & filename, const Mesh & mesh) {
std::ofstream outfile;
const Vector<Real> & nodes = mesh.getNodes();
outfile.open(filename.c_str());
outfile << "$MeshFormat" << std::endl;
outfile << "2 0 8" << std::endl;;
outfile << "$EndMeshFormat" << std::endl;;
outfile << "$Nodes" << std::endl;;
outfile << nodes.getSize() << std::endl;;
outfile << std::uppercase;
for(UInt i = 0; i < nodes.getSize(); ++i) {
Int offset = i * nodes.getNbComponent();
outfile << i+1;
for(UInt j = 0; j < nodes.getNbComponent(); ++j) {
outfile << " " << nodes.values[offset + j];
}
if(nodes.getNbComponent() == 2)
outfile << " " << 0.;
outfile << std::endl;;
}
outfile << std::nouppercase;
outfile << "$EndNodes" << std::endl;;
outfile << "$Elements" << std::endl;;
const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList();
Mesh::ConnectivityTypeList::const_iterator it;
Int nb_elements = 0;
for(it = type_list.begin(); it != type_list.end(); ++it) {
const Vector<UInt> & connectivity = mesh.getConnectivity(*it);
nb_elements += connectivity.getSize();
}
outfile << nb_elements << std::endl;
UInt element_idx = 1;
for(it = type_list.begin(); it != type_list.end(); ++it) {
ElementType type = *it;
const Vector<UInt> & connectivity = mesh.getConnectivity(type);
for(UInt i = 0; i < connectivity.getSize(); ++i) {
UInt offset = i * connectivity.getNbComponent();
outfile << element_idx << " " << type << " 3 0 0 0";
for(UInt j = 0; j < connectivity.getNbComponent(); ++j) {
outfile << " " << connectivity.values[offset + j];
}
outfile << std::endl;
}
element_idx++;
}
outfile << "$EndElements" << std::endl;;
outfile.close();
}
/* -------------------------------------------------------------------------- */
__END_AKANTU__
diff --git a/mesh_utils/mesh_partition.hh b/mesh_utils/mesh_partition.hh
index 7cc3e001b..4d89c59ef 100644
--- a/mesh_utils/mesh_partition.hh
+++ b/mesh_utils/mesh_partition.hh
@@ -1,109 +1,109 @@
/**
- * @file mesh_partition.h
+ * @file mesh_partition.hh
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Thu Aug 12 16:24:40 2010
*
* @brief tools to partitionate a mesh
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -------------------------------------------------------------------------- */
-#ifndef __AKANTU_MESH_PARTITION_H__
-#define __AKANTU_MESH_PARTITION_H__
+#ifndef __AKANTU_MESH_PARTITION_HH__
+#define __AKANTU_MESH_PARTITION_HH__
/* -------------------------------------------------------------------------- */
#include "aka_memory.hh"
#include "mesh.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
class MeshPartition : public Memory {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
MeshPartition(const Mesh & mesh, UInt spatial_dimension,
const MemoryID & memory_id = 0);
virtual ~MeshPartition();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
virtual void partitionate(UInt nb_part) = 0;
/// function to print the contain of the class
//virtual void printself(std::ostream & stream, int indent = 0) const = 0;
protected:
/// build the dual graph of the mesh, for all element of spatial_dimension
void buildDualGraph(Vector<Int> & dxadj, Vector<Int> & dadjncy);
/// fill the partitions array with a given linearized partition information
void fillPartitionInformations(const Mesh & mesh, const Int * linearized_partitions);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Partition, partitions, const Vector<UInt> &);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE(GhostPartition,
ghost_partitions, const Vector<UInt> &);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE(GhostPartitionOffset,
ghost_partitions_offset, const Vector<UInt> &);
AKANTU_GET_MACRO(NbPartition, nb_partitions, UInt);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// id
std::string id;
/// the mesh to partition
const Mesh & mesh;
/// dimension of the elements to consider in the mesh
UInt spatial_dimension;
/// number of partitions
UInt nb_partitions;
/// partition numbers
ByElementTypeUInt partitions;
ByElementTypeUInt ghost_partitions;
ByElementTypeUInt ghost_partitions_offset;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
//#include "mesh_partition_inline_impl.cc"
/// standard output stream operator
// inline std::ostream & operator <<(std::ostream & stream, const MeshPartition & _this)
// {
// _this.printself(stream);
// return stream;
// }
__END_AKANTU__
-#endif /* __AKANTU_MESH_PARTITION_H__ */
+#endif /* __AKANTU_MESH_PARTITION_HH__ */
diff --git a/mesh_utils/mesh_partition/mesh_partition_scotch.cc b/mesh_utils/mesh_partition/mesh_partition_scotch.cc
index 02a12639f..0ea081bf1 100644
--- a/mesh_utils/mesh_partition/mesh_partition_scotch.cc
+++ b/mesh_utils/mesh_partition/mesh_partition_scotch.cc
@@ -1,169 +1,169 @@
/**
* @file mesh_partition_scotch.cc
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Fri Aug 13 11:54:11 2010
*
* @brief implementation of the MeshPartitionScotch class
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -------------------------------------------------------------------------- */
#include <cstdio>
#include <fstream>
-extern "C" {
+//extern "C" {
#include <scotch.h>
-}
+//}
/* -------------------------------------------------------------------------- */
#include "mesh_partition_scotch.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
MeshPartitionScotch::MeshPartitionScotch(const Mesh & mesh, UInt spatial_dimension,
const MemoryID & memory_id) :
MeshPartition(mesh, spatial_dimension, memory_id) {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MeshPartitionScotch::partitionate(UInt nb_part) {
AKANTU_DEBUG_IN();
nb_partitions = nb_part;
AKANTU_DEBUG_INFO("Partitioning the mesh " << mesh.getID()
<< " in " << nb_part << " parts.");
Vector<Int> dxadj;
Vector<Int> dadjncy;
buildDualGraph(dxadj, dadjncy);
/// variables that will hold our structures in scotch format
SCOTCH_Graph scotch_graph;
SCOTCH_Strat scotch_strat;
/// description number and arrays for struct mesh for scotch
SCOTCH_Num baseval = 0; //base numbering for element and
//nodes (0 -> C , 1 -> fortran)
SCOTCH_Num vertnbr = dxadj.getSize() - 1; //number of vertexes
SCOTCH_Num *parttab; //array of partitions
SCOTCH_Num edgenbr = dxadj.values[vertnbr]; //twice the number of "edges"
//(an "edge" bounds two nodes)
SCOTCH_Num * verttab = dxadj.values; //array of start indices in edgetab
SCOTCH_Num * vendtab = NULL; //array of after-last indices in edgetab
SCOTCH_Num * velotab = NULL; //integer load associated with
//every vertex ( optional )
- SCOTCH_Num *edlotab = NULL; //integer load associated with
+ SCOTCH_Num * edlotab = NULL; //integer load associated with
//every edge ( optional )
- SCOTCH_Num *edgetab = dadjncy.values; // adjacency array of every vertex
- SCOTCH_Num *vlbltab = NULL; // vertex label array (optional)
+ SCOTCH_Num * edgetab = dadjncy.values; // adjacency array of every vertex
+ SCOTCH_Num * vlbltab = NULL; // vertex label array (optional)
/// Allocate space for Scotch arrays
parttab = new SCOTCH_Num[vertnbr];
/// Initialize the strategy structure
SCOTCH_stratInit (&scotch_strat);
/// Initialize the graph structure
SCOTCH_graphInit(&scotch_graph);
/// Build the graph from the adjacency arrays
SCOTCH_graphBuild(&scotch_graph, baseval,
vertnbr, verttab, vendtab,
velotab, vlbltab, edgenbr,
edgetab, edlotab);
/// Check the graph
AKANTU_DEBUG_ASSERT(SCOTCH_graphCheck(&scotch_graph) == 0,
"Graph to partition is not consistent");
#ifndef AKANTU_NDEBUG
if (AKANTU_DEBUG_TEST(dblDump)) {
/// save initial graph
FILE *fgraphinit = fopen("GraphIniFile.grf", "w");
SCOTCH_graphSave(&scotch_graph,fgraphinit);
fclose(fgraphinit);
/// write geometry file
std::ofstream fgeominit;
fgeominit.open("GeomIniFile.xyz");
fgeominit << spatial_dimension << std::endl << vertnbr << std::endl;
const Vector<Real> & nodes = mesh.getNodes();
const Mesh::ConnectivityTypeList & f_type_list = mesh.getConnectivityTypeList();
Mesh::ConnectivityTypeList::const_iterator f_it;
UInt out_linerized_el = 0;
for(f_it = f_type_list.begin(); f_it != f_type_list.end(); ++f_it) {
ElementType type = *f_it;
if(Mesh::getSpatialDimension(type) != mesh.getSpatialDimension()) continue;
UInt nb_element = mesh.getNbElement(*f_it);
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
const Vector<UInt> & connectivity = mesh.getConnectivity(type);
Real mid[spatial_dimension] ;
for (UInt el = 0; el < nb_element; ++el) {
memset(mid, 0, spatial_dimension*sizeof(Real));
for (UInt n = 0; n < nb_nodes_per_element; ++n) {
UInt node = connectivity.values[nb_nodes_per_element * el + n];
for (UInt s = 0; s < spatial_dimension; ++s)
mid[s] += ((Real) nodes.values[node * spatial_dimension + s]) / ((Real) nb_nodes_per_element);
}
fgeominit << out_linerized_el++ << " ";
for (UInt s = 0; s < spatial_dimension; ++s)
fgeominit << mid[s] << " ";
fgeominit << std::endl;;
}
}
fgeominit.close();
}
#endif
/// Partition the mesh
SCOTCH_graphPart(&scotch_graph, nb_part, &scotch_strat, parttab);
/// Check the graph
AKANTU_DEBUG_ASSERT(SCOTCH_graphCheck(&scotch_graph) == 0,
"Partitioned graph is not consistent");
#ifndef AKANTU_NDEBUG
if (AKANTU_DEBUG_TEST(dblDump)) {
/// save the partitioned graph
FILE *fgraph=fopen("GraphFile.grf", "w");
SCOTCH_graphSave(&scotch_graph, fgraph);
fclose(fgraph);
/// save the partition map
std::ofstream fmap;
fmap.open("MapFile.map");
fmap << vertnbr << std::endl;
for (Int i = 0; i < vertnbr; i++) fmap << i << " " << parttab[i] << std::endl;
fmap.close();
}
#endif
/// free the scotch data structures
SCOTCH_stratExit(&scotch_strat);
SCOTCH_graphFree(&scotch_graph);
SCOTCH_graphExit(&scotch_graph);
fillPartitionInformations(mesh, parttab);
delete [] parttab;
AKANTU_DEBUG_OUT();
}
__END_AKANTU__
diff --git a/mesh_utils/mesh_utils.cc b/mesh_utils/mesh_utils.cc
index b6b75d272..d8336bda0 100644
--- a/mesh_utils/mesh_utils.cc
+++ b/mesh_utils/mesh_utils.cc
@@ -1,475 +1,475 @@
/**
* @file mesh_utils.cc
* @author Guillaume ANCIAUX <anciaux@lsmscluster1.epfl.ch>
* @date Wed Aug 18 14:21:00 2010
*
* @brief All mesh utils necessary for various tasks
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -------------------------------------------------------------------------- */
#include "mesh_utils.hh"
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
void MeshUtils::buildNode2Elements(const Mesh & mesh,
Vector<UInt> & node_offset,
Vector<UInt> & node_to_elem,
UInt spatial_dimension) {
AKANTU_DEBUG_IN();
if (spatial_dimension == 0) spatial_dimension = mesh.getSpatialDimension();
UInt nb_nodes = mesh.getNbNodes();
const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList();
Mesh::ConnectivityTypeList::const_iterator it;
UInt nb_types = type_list.size();
UInt nb_good_types = 0;
UInt nb_nodes_per_element[nb_types];
UInt nb_nodes_per_element_p1[nb_types];
UInt * conn_val[nb_types];
UInt nb_element[nb_types];
ElementType type_p1;
for(it = type_list.begin(); it != type_list.end(); ++it) {
ElementType type = *it;
if(Mesh::getSpatialDimension(type) != spatial_dimension) continue;
nb_nodes_per_element[nb_good_types] = Mesh::getNbNodesPerElement(type);
type_p1 = Mesh::getP1ElementType(type);
nb_nodes_per_element_p1[nb_good_types] = Mesh::getNbNodesPerElement(type_p1);
conn_val[nb_good_types] = mesh.getConnectivity(type).values;
nb_element[nb_good_types] = mesh.getConnectivity(type).getSize();
nb_good_types++;
}
AKANTU_DEBUG_ASSERT(nb_good_types != 0,
"Some elements must be found in right dimension to compute facets!");
/// array for the node-element list
node_offset.resize(nb_nodes + 1);
UInt * node_offset_val = node_offset.values;
/// count number of occurrence of each node
memset(node_offset_val, 0, (nb_nodes + 1)*sizeof(UInt));
for (UInt t = 0; t < nb_good_types; ++t) {
for (UInt el = 0; el < nb_element[t]; ++el) {
UInt el_offset = el*nb_nodes_per_element[t];
for (UInt n = 0; n < nb_nodes_per_element_p1[t]; ++n) {
node_offset_val[conn_val[t][el_offset + n]]++;
}
}
}
/// convert the occurrence array in a csr one
for (UInt i = 1; i < nb_nodes; ++i) node_offset_val[i] += node_offset_val[i-1];
for (UInt i = nb_nodes; i > 0; --i) node_offset_val[i] = node_offset_val[i-1];
node_offset_val[0] = 0;
/// rearrange element to get the node-element list
node_to_elem.resize(node_offset_val[nb_nodes]);
UInt * node_to_elem_val = node_to_elem.values;
for (UInt t = 0, linearized_el = 0; t < nb_good_types; ++t)
for (UInt el = 0; el < nb_element[t]; ++el, ++linearized_el) {
UInt el_offset = el*nb_nodes_per_element[t];
for (UInt n = 0; n < nb_nodes_per_element_p1[t]; ++n)
node_to_elem_val[node_offset_val[conn_val[t][el_offset + n]]++] = linearized_el;
}
for (UInt i = nb_nodes; i > 0; --i) node_offset_val[i] = node_offset_val[i-1];
node_offset_val[0] = 0;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MeshUtils::buildNode2ElementsByElementType(const Mesh & mesh,
ElementType type,
Vector<UInt> & node_offset,
Vector<UInt> & node_to_elem) {
AKANTU_DEBUG_IN();
UInt nb_nodes = mesh.getNbNodes();
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
UInt nb_elements = mesh.getConnectivity(type).getSize();
UInt * conn_val = mesh.getConnectivity(type).values;
/// array for the node-element list
node_offset.resize(nb_nodes + 1);
UInt * node_offset_val = node_offset.values;
memset(node_offset_val, 0, (nb_nodes + 1)*sizeof(UInt));
/// count number of occurrence of each node
for (UInt el = 0; el < nb_elements; ++el)
for (UInt n = 0; n < nb_nodes_per_element; ++n)
node_offset_val[conn_val[nb_nodes_per_element*el + n]]++;
/// convert the occurrence array in a csr one
for (UInt i = 1; i < nb_nodes; ++i) node_offset_val[i] += node_offset_val[i-1];
for (UInt i = nb_nodes; i > 0; --i) node_offset_val[i] = node_offset_val[i-1];
node_offset_val[0] = 0;
/// save the element index in the node-element list
node_to_elem.resize(node_offset_val[nb_nodes]);
UInt * node_to_elem_val = node_to_elem.values;
for (UInt el = 0; el < nb_elements; ++el)
for (UInt n = 0; n < nb_nodes_per_element; ++n) {
node_to_elem_val[node_offset_val[conn_val[nb_nodes_per_element*el + n]]] = el;
node_offset_val[conn_val[nb_nodes_per_element*el + n]]++;
}
/// rearrange node_offset to start with 0
for (UInt i = nb_nodes; i > 0; --i) node_offset_val[i] = node_offset_val[i-1];
node_offset_val[0] = 0;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MeshUtils::buildFacets(Mesh & mesh, bool boundary_flag, bool internal_flag){
AKANTU_DEBUG_IN();
Vector<UInt> node_offset;
Vector<UInt> node_to_elem;
buildNode2Elements(mesh,node_offset,node_to_elem);
// std::cout << "node offset " << std::endl << node_offset << std::endl;
// std::cout << "node to elem " << std::endl << node_to_elem << std::endl;
const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList();
Mesh::ConnectivityTypeList::const_iterator it;
UInt nb_types = type_list.size();
UInt nb_good_types = 0;
UInt nb_nodes_per_element[nb_types];
UInt nb_nodes_per_facet[nb_types];
UInt nb_facets[nb_types];
UInt ** node_in_facet[nb_types];
Vector<UInt> * connectivity_facets[nb_types];
Vector<UInt> * connectivity_internal_facets[nb_types];
UInt * conn_val[nb_types];
UInt nb_element[nb_types];
ElementType facet_type;
Mesh * internal_facets_mesh = NULL;
UInt spatial_dimension = mesh.getSpatialDimension();
if (internal_flag) internal_facets_mesh = mesh.getInternalFacetsMeshPointer();
for(it = type_list.begin(); it != type_list.end(); ++it) {
ElementType type = *it;
if(mesh.getSpatialDimension(type) != spatial_dimension) continue;
nb_nodes_per_element[nb_good_types] = mesh.getNbNodesPerElement(type);
facet_type = mesh.getFacetElementType(type);
nb_facets[nb_good_types] = mesh.getNbFacetsPerElement(type);
node_in_facet[nb_good_types] = mesh.getFacetLocalConnectivity(type);
nb_nodes_per_facet[nb_good_types] = mesh.getNbNodesPerElement(facet_type);
if (boundary_flag){
// getting connectivity of boundary facets
connectivity_facets[nb_good_types] = mesh.getConnectivityPointer(facet_type);
connectivity_facets[nb_good_types]->resize(0);
}
if (internal_flag){
// getting connectivity of internal facets
connectivity_internal_facets[nb_good_types] =
internal_facets_mesh->getConnectivityPointer(facet_type);
connectivity_internal_facets[nb_good_types]->resize(0);
}
conn_val[nb_good_types] = mesh.getConnectivity(type).values;
nb_element[nb_good_types] = mesh.getConnectivity(type).getSize();
nb_good_types++;
}
Vector<UInt> counter;
/// count number of occurrence of each node
for (UInt t = 0,linearized_el = 0; t < nb_good_types; ++t) {
for (UInt el = 0; el < nb_element[t]; ++el, ++linearized_el) {
UInt el_offset = el*nb_nodes_per_element[t];
for (UInt f = 0; f < nb_facets[t]; ++f) {
//build the nodes involved in facet 'f'
UInt facet_nodes[nb_nodes_per_facet[t]];
for (UInt n = 0; n < nb_nodes_per_facet[t]; ++n) {
UInt node_facet = node_in_facet[t][f][n];
facet_nodes[n] = conn_val[t][el_offset + node_facet];
}
//our reference is the first node
UInt * first_node_elements = node_to_elem.values+node_offset.values[facet_nodes[0]];
UInt first_node_nelements =
node_offset.values[facet_nodes[0]+1]-
node_offset.values[facet_nodes[0]];
counter.resize(first_node_nelements);
memset(counter.values,0,sizeof(UInt)*first_node_nelements);
//loop over the other nodes to search intersecting elements
for (UInt n = 1; n < nb_nodes_per_facet[t]; ++n) {
UInt * node_elements = node_to_elem.values+node_offset.values[facet_nodes[n]];
UInt node_nelements =
node_offset.values[facet_nodes[n]+1]-
node_offset.values[facet_nodes[n]];
for (UInt el1 = 0; el1 < first_node_nelements; ++el1) {
for (UInt el2 = 0; el2 < node_nelements; ++el2) {
if (first_node_elements[el1] == node_elements[el2]) {
++counter.values[el1];
break;
}
}
if (counter.values[el1] == nb_nodes_per_facet[t]) break;
}
}
bool connected_facet = false;
//the connected elements are those for which counter is n_facet
// UInt connected_element = -1;
for (UInt el1 = 0; el1 < counter.getSize(); el1++) {
UInt el_index = node_to_elem.values[node_offset.values[facet_nodes[0]]+el1];
if (counter.values[el1] == nb_nodes_per_facet[t]-1 && el_index > linearized_el){
// connected_element = el_index;
AKANTU_DEBUG(dblDump,"connecting elements " << linearized_el << " and " << el_index);
if (internal_flag)
connectivity_internal_facets[t]->push_back(facet_nodes);
}
if (counter.values[el1] == nb_nodes_per_facet[t]-1 && el_index != linearized_el)
connected_facet = true;
}
if (!connected_facet) {
AKANTU_DEBUG(dblDump,"facet " << f << " in element " << linearized_el << " is a boundary");
if (boundary_flag)
connectivity_facets[t]->push_back(facet_nodes);
}
}
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
-void MeshUtils::buildNormals(Mesh & mesh,UInt spatial_dimension){
- AKANTU_DEBUG_IN();
- const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList();
- Mesh::ConnectivityTypeList::const_iterator it;
-
- UInt nb_types = type_list.size();
- UInt nb_nodes_per_element[nb_types];
- UInt nb_nodes_per_element_p1[nb_types];
-
- UInt nb_good_types = 0;
-
- Vector<UInt> * connectivity[nb_types];
- Vector<Real> * normals[nb_types];
-
- if (spatial_dimension == 0) spatial_dimension = mesh.getSpatialDimension();
-
- for(it = type_list.begin(); it != type_list.end(); ++it) {
- ElementType type = *it;
- ElementType type_p1 = mesh.getP1ElementType(type);
- if(mesh.getSpatialDimension(type) != spatial_dimension) continue;
-
- nb_nodes_per_element[nb_good_types] = mesh.getNbNodesPerElement(type);
- nb_nodes_per_element_p1[nb_good_types] = mesh.getNbNodesPerElement(type_p1);
-
- // getting connectivity
- connectivity[nb_good_types] = mesh.getConnectivityPointer(type);
- if (!connectivity[nb_good_types])
- AKANTU_DEBUG_ERROR("connectivity is not allocatted : this should probably not have happened");
-
- //getting array of normals
- normals[nb_good_types] = mesh.getNormalsPointer(type);
- if(normals[nb_good_types])
- normals[nb_good_types]->resize(0);
- else
- normals[nb_good_types] = &mesh.createNormals(type);
-
- nb_good_types++;
- }
- AKANTU_DEBUG_OUT();
-}
+// void MeshUtils::buildNormals(Mesh & mesh,UInt spatial_dimension){
+// AKANTU_DEBUG_IN();
+// const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList();
+// Mesh::ConnectivityTypeList::const_iterator it;
+
+// UInt nb_types = type_list.size();
+// UInt nb_nodes_per_element[nb_types];
+// UInt nb_nodes_per_element_p1[nb_types];
+
+// UInt nb_good_types = 0;
+
+// Vector<UInt> * connectivity[nb_types];
+// Vector<Real> * normals[nb_types];
+
+// if (spatial_dimension == 0) spatial_dimension = mesh.getSpatialDimension();
+
+// for(it = type_list.begin(); it != type_list.end(); ++it) {
+// ElementType type = *it;
+// ElementType type_p1 = mesh.getP1ElementType(type);
+// if(mesh.getSpatialDimension(type) != spatial_dimension) continue;
+
+// nb_nodes_per_element[nb_good_types] = mesh.getNbNodesPerElement(type);
+// nb_nodes_per_element_p1[nb_good_types] = mesh.getNbNodesPerElement(type_p1);
+
+// // getting connectivity
+// connectivity[nb_good_types] = mesh.getConnectivityPointer(type);
+// if (!connectivity[nb_good_types])
+// AKANTU_DEBUG_ERROR("connectivity is not allocatted : this should probably not have happened");
+
+// //getting array of normals
+// normals[nb_good_types] = mesh.getNormalsPointer(type);
+// if(normals[nb_good_types])
+// normals[nb_good_types]->resize(0);
+// else
+// normals[nb_good_types] = &mesh.createNormals(type);
+
+// nb_good_types++;
+// }
+// AKANTU_DEBUG_OUT();
+// }
/* -------------------------------------------------------------------------- */
void MeshUtils::renumberMeshNodes(Mesh & mesh,
UInt * local_connectivities,
UInt nb_local_element,
UInt nb_ghost_element,
ElementType type,
Vector<UInt> & nodes_numbers) {
AKANTU_DEBUG_IN();
nodes_numbers.resize(0);
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
UInt * conn = local_connectivities;
/// renumber the nodes
for (UInt el = 0; el < nb_local_element + nb_ghost_element; ++el) {
for (UInt n = 0; n < nb_nodes_per_element; ++n) {
Int nn = nodes_numbers.find(*conn);
if(nn == -1) {
nodes_numbers.push_back(*conn);
*conn = nodes_numbers.getSize() - 1;
} else {
*conn = nn;
}
*conn++;
}
}
/// copy the renumbered connectivity to the right place
Vector<UInt> * local_conn = mesh.getConnectivityPointer(type);
local_conn->resize(nb_local_element);
memcpy(local_conn->values,
local_connectivities,
nb_local_element * nb_nodes_per_element * sizeof(UInt));
Vector<UInt> * ghost_conn = mesh.getGhostConnectivityPointer(type);
ghost_conn->resize(nb_ghost_element);
memcpy(ghost_conn->values,
local_connectivities + nb_local_element * nb_nodes_per_element,
nb_ghost_element * nb_nodes_per_element * sizeof(UInt));
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void MeshUtils::buildSurfaceID(Mesh & mesh) {
AKANTU_DEBUG_IN();
Vector<UInt> node_offset;
Vector<UInt> node_to_elem;
/// Get list of surface elements
UInt spatial_dimension = mesh.getSpatialDimension();
buildNode2Elements(mesh, node_offset, node_to_elem, spatial_dimension-1);
/// Find which types of elements have been linearized
const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList();
Mesh::ConnectivityTypeList::const_iterator it;
UInt nb_types = type_list.size();
ElementType lin_element_type[nb_types];
UInt nb_lin_types = 0;
UInt nb_nodes_per_element[nb_types];
UInt nb_nodes_per_element_p1[nb_types];
UInt * conn_val[nb_types];
UInt nb_element[nb_types];
ElementType type_p1;
for(it = type_list.begin(); it != type_list.end(); ++it) {
ElementType type = *it;
if(Mesh::getSpatialDimension(type) != spatial_dimension) continue;
ElementType facet_type = mesh.getFacetElementType(type);
lin_element_type[nb_lin_types] = facet_type;
nb_nodes_per_element[nb_lin_types] = Mesh::getNbNodesPerElement(facet_type);
type_p1 = Mesh::getP1ElementType(facet_type);
nb_nodes_per_element_p1[nb_lin_types] = Mesh::getNbNodesPerElement(type_p1);
conn_val[nb_lin_types] = mesh.getConnectivity(facet_type).values;
nb_element[nb_lin_types] = mesh.getConnectivity(facet_type).getSize();
nb_lin_types++;
}
for (UInt i = 1; i < nb_lin_types; ++i) nb_element[i] += nb_element[i+1];
for (UInt i = nb_lin_types; i > 0; --i) nb_element[i] = nb_element[i-1];
nb_element[0] = 0;
/// Find close surfaces
Vector<Int> surface_value_id(1, nb_element[nb_lin_types], -1);
Int * surf_val = surface_value_id.values;
UInt nb_surfaces = 0;
UInt nb_cecked_elements;
UInt nb_elements_to_ceck;
UInt * elements_to_ceck = new UInt [nb_element[nb_lin_types]];
memset(elements_to_ceck, 0, nb_element[nb_lin_types]*sizeof(UInt));
for (UInt lin_el = 0; lin_el < nb_element[nb_lin_types]; ++lin_el) {
if(surf_val[lin_el] != -1) continue; /* Surface id already assigned */
/* First element of new surface */
surf_val[lin_el] = nb_surfaces;
nb_cecked_elements = 0;
nb_elements_to_ceck = 1;
memset(elements_to_ceck, 0, nb_element[nb_lin_types]*sizeof(UInt));
elements_to_ceck[0] = lin_el;
// Find others elements belonging to this surface
while(nb_cecked_elements < nb_elements_to_ceck) {
UInt ceck_lin_el = elements_to_ceck[nb_cecked_elements];
// Transform linearized index of element into ElementType one
UInt lin_type_nb = 0;
while (ceck_lin_el >= nb_element[lin_type_nb+1])
lin_type_nb++;
UInt ceck_el = ceck_lin_el - nb_element[lin_type_nb];
// Get connected elements
UInt el_offset = ceck_el*nb_nodes_per_element[lin_type_nb];
for (UInt n = 0; n < nb_nodes_per_element_p1[lin_type_nb]; ++n) {
UInt node_id = conn_val[lin_type_nb][el_offset + n];
for (UInt i = node_offset.values[node_id]; i < node_offset.values[node_id+1]; ++i)
if(surf_val[node_to_elem.values[i]] == -1) { /* Found new surface element */
surf_val[node_to_elem.values[i]] = nb_surfaces;
elements_to_ceck[nb_elements_to_ceck] = node_to_elem.values[i];
nb_elements_to_ceck++;
}
}
nb_cecked_elements++;
}
nb_surfaces++;
}
delete [] elements_to_ceck;
/// Transform local linearized element index in the global one
for (UInt i = 0; i < nb_lin_types; ++i) nb_element[i] = nb_element[i+1] - nb_element[i];
UInt el_offset = 0;
for (UInt type_it = 0; type_it < nb_lin_types; ++type_it) {
ElementType type = lin_element_type[type_it];
Vector<UInt> * surf_id_type = mesh.getSurfaceIdPointer(type);
surf_id_type->resize(nb_element[type_it]);
for (UInt el = 0; el < nb_element[type_it]; ++el)
surf_id_type->values[el] = surf_val[el+el_offset];
el_offset += nb_element[type_it];
}
/// Set nb_surfaces in mesh
mesh.nb_surfaces = nb_surfaces;
AKANTU_DEBUG_OUT();
}
__END_AKANTU__
// LocalWords: ElementType
diff --git a/mesh_utils/mesh_utils.hh b/mesh_utils/mesh_utils.hh
index 567e259d5..721112367 100644
--- a/mesh_utils/mesh_utils.hh
+++ b/mesh_utils/mesh_utils.hh
@@ -1,91 +1,91 @@
/**
* @file mesh_utils.hh
* @author Guillaume ANCIAUX <anciaux@lsmscluster1.epfl.ch>
* @date Wed Aug 18 14:03:39 2010
*
* @brief All mesh utils necessary for various tasks
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MESH_UTILS_HH__
#define __AKANTU_MESH_UTILS_HH__
#include "aka_common.hh"
#include "mesh.hh"
__BEGIN_AKANTU__
class MeshUtils {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
MeshUtils();
virtual ~MeshUtils();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// build map from nodes to elements
static void buildNode2Elements(const Mesh & mesh, Vector<UInt> & node_offset, Vector<UInt> & node_to_elem, UInt spatial_dimension = 0);
/// build map from nodes to elements for a specific element type
static void buildNode2ElementsByElementType(const Mesh & mesh, ElementType type, Vector<UInt> & node_offset, Vector<UInt> & node_to_elem);
/// build facets elements : boundary and/or internals
static void buildFacets(Mesh & mesh, bool boundary_flag=1, bool internal_flag=0);
/// build normal to some elements
- static void buildNormals(Mesh & mesh, UInt spatial_dimension=0);
+ // static void buildNormals(Mesh & mesh, UInt spatial_dimension=0);
/// take the local_connectivity array as the array of local and ghost
/// connectivity, renumber the nodes and set the connectivity of the mesh
static void renumberMeshNodes(Mesh & mesh,
UInt * local_connectivities,
UInt nb_local_element,
UInt nb_ghost_element,
ElementType type,
Vector<UInt> & old_nodes);
/// Detect closed surfaces of the mesh and save the surface id
/// of the surface elements in the array surface_id
static void buildSurfaceID(Mesh & mesh);
/// function to print the contain of the class
// virtual void printself(std::ostream & stream, int indent = 0) const;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#include "mesh_utils_inline_impl.cc"
/// standard output stream operator
// inline std::ostream & operator <<(std::ostream & stream, const MeshUtils & _this)
// {
// _this.printself(stream);
// return stream;
// }
__END_AKANTU__
#endif /* __AKANTU_MESH_UTILS_HH__ */
diff --git a/model/model_inline_impl.cc b/model/model_inline_impl.cc
index 79b5ada89..bae51cea7 100644
--- a/model/model_inline_impl.cc
+++ b/model/model_inline_impl.cc
@@ -1,45 +1,50 @@
/**
* @file model_inline_impl.cc
* @author Guillaume ANCIAUX <anciaux@lsmscluster1.epfl.ch>
* @date Fri Aug 20 17:18:08 2010
*
* @brief inline implementation of the model class
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -------------------------------------------------------------------------- */
inline Model::Model(Mesh & mesh,
UInt spatial_dimension,
const ModelID & id,
const MemoryID & memory_id) :
Memory(memory_id), id(id) {
AKANTU_DEBUG_IN();
this->spatial_dimension = (spatial_dimension == 0) ? mesh.getSpatialDimension() : spatial_dimension;
std::stringstream sstr; sstr << id << ":fem";
this->fem = new FEM(mesh, spatial_dimension, sstr.str(), memory_id);
this->fem_boundary = NULL;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
inline Model::~Model() {
AKANTU_DEBUG_IN();
delete fem;
+
+ if(fem_boundary) delete fem_boundary;
+
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
inline FEM & Model::getFEMBoundary(){
AKANTU_DEBUG_IN();
+
if (!fem_boundary){
MeshUtils::buildFacets(fem->getMesh());
- std::stringstream sstr; sstr << id << ":femboundary";
- this->fem_boundary = new FEM(fem->getMesh(), spatial_dimension-1, sstr.str(), memory_id);
+ std::stringstream sstr; sstr << id << ":femboundary";
+ this->fem_boundary = new FEM(fem->getMesh(), spatial_dimension-1, sstr.str(), memory_id);
}
+
AKANTU_DEBUG_OUT();
return *this->fem_boundary;
}
diff --git a/model/solid_mechanics/material.cc b/model/solid_mechanics/material.cc
index aad6bef71..d2bd3ede3 100644
--- a/model/solid_mechanics/material.cc
+++ b/model/solid_mechanics/material.cc
@@ -1,252 +1,252 @@
/**
* @file material.cc
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Tue Jul 27 11:43:41 2010
*
* @brief Implementation of the common part of the material class
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -------------------------------------------------------------------------- */
#include "material.hh"
#include "solid_mechanics_model.hh"
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
Material::Material(SolidMechanicsModel & model, const MaterialID & id) :
Memory(model.getMemoryID()), id(id),
spatial_dimension(model.getSpatialDimension()), name(""),
model(&model), potential_energy_vector(false),
is_init(false) {
AKANTU_DEBUG_IN();
for(UInt t = _not_defined; t < _max_element_type; ++t) {
this->stress [t] = NULL;
this->strain [t] = NULL;
this->potential_energy [t] = NULL;
this->element_filter [t] = NULL;
this->ghost_stress [t] = NULL;
this->ghost_strain [t] = NULL;
this->ghost_potential_energy[t] = NULL;
this->ghost_element_filter [t] = NULL;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
Material::~Material() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void Material::setParam(const std::string & key, const std::string & value,
__attribute__ ((unused)) const MaterialID & id) {
if(key == "name") name = std::string(value);
else AKANTU_DEBUG_ERROR("Unknown material property : " << key);
}
/* -------------------------------------------------------------------------- */
void Material::initMaterial() {
AKANTU_DEBUG_IN();
/// for each connectivity types allocate the element filer array of the material
UInt spatial_dimension = model->getSpatialDimension();
const Mesh::ConnectivityTypeList & type_list =
model->getFEM().getMesh().getConnectivityTypeList();
Mesh::ConnectivityTypeList::const_iterator it;
for(it = type_list.begin(); it != type_list.end(); ++it) {
if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue;
std::stringstream sstr; sstr << id << ":element_filer:"<< *it;
element_filter[*it] = &(alloc<UInt> (sstr.str(), 0, 1));
std::stringstream sstr_stre; sstr_stre << id << ":stress:" << *it;
std::stringstream sstr_stra; sstr_stra << id << ":strain:" << *it;
stress[*it] = &(alloc<Real>(sstr_stre.str(), 0,
spatial_dimension * spatial_dimension));
strain[*it] = &(alloc<Real>(sstr_stra.str(), 0,
spatial_dimension * spatial_dimension));
}
const Mesh::ConnectivityTypeList & ghost_type_list =
model->getFEM().getMesh().getConnectivityTypeList(_ghost);
for(it = ghost_type_list.begin(); it != ghost_type_list.end(); ++it) {
if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue;
std::stringstream sstr; sstr << id << ":ghost_element_filer:"<< *it;
ghost_element_filter[*it] = &(alloc<UInt> (sstr.str(), 0, 1));
std::stringstream sstr_stre; sstr_stre << id << ":ghost_stress:" << *it;
std::stringstream sstr_stra; sstr_stra << id << ":ghost_strain:" << *it;
ghost_stress[*it] = &(alloc<Real>(sstr_stre.str(), 0,
spatial_dimension * spatial_dimension));
ghost_strain[*it] = &(alloc<Real>(sstr_stra.str(), 0,
spatial_dimension * spatial_dimension));
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Compute the residual by assembling @f$\int_{e} \sigma_e \frac{\partial
* \varphi}{\partial X} dX @f$
*
* @param[in] current_position nodes postition + displacements
* @param[in] ghost_type compute the residual for _ghost or _not_ghost element
*/
void Material::updateResidual(Vector<Real> & current_position, GhostType ghost_type) {
AKANTU_DEBUG_IN();
UInt spatial_dimension = model->getSpatialDimension();
Vector<Real> & residual = const_cast<Vector<Real> &>(model->getResidual());
const Mesh::ConnectivityTypeList & type_list =
model->getFEM().getMesh().getConnectivityTypeList(ghost_type);
Mesh::ConnectivityTypeList::const_iterator it;
for(it = type_list.begin(); it != type_list.end(); ++it) {
if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue;
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it);
UInt size_of_shapes_derivatives = FEM::getShapeDerivativesSize(*it);
UInt nb_quadrature_points = FEM::getNbQuadraturePoints(*it);
Vector<Real> * strain_vect;
Vector<Real> * stress_vect;
Vector<UInt> * elem_filter;
const Vector<Real> * shapes_derivatives;
if(ghost_type == _not_ghost) {
elem_filter = element_filter[*it];
strain_vect = strain[*it];
stress_vect = stress[*it];
shapes_derivatives = &(model->getFEM().getShapesDerivatives(*it));
} else {
elem_filter = ghost_element_filter[*it];
strain_vect = ghost_strain[*it];
stress_vect = ghost_stress[*it];
shapes_derivatives = &(model->getFEM().getGhostShapesDerivatives(*it));
}
UInt nb_element = elem_filter->getSize();
/// compute @f$\nabla u@f$
model->getFEM().gradientOnQuadraturePoints(current_position, *strain_vect,
spatial_dimension,
*it, ghost_type, elem_filter);
/// compute @f$\mathbf{\sigma}_q@f$ from @f$\nabla u@f$
computeStress(*it, ghost_type);
/// compute @f$\sigma \frac{\partial \varphi}{\partial X}@f$ by @f$\mathbf{B}^t \mathbf{\sigma}_q@f$
Vector<Real> * sigma_dphi_dx =
new Vector<Real>(nb_element*nb_quadrature_points, size_of_shapes_derivatives, "sigma_x_dphi_/_dX");
Real * shapesd = shapes_derivatives->values;
UInt size_of_shapesd = shapes_derivatives->getNbComponent();
Real * shapesd_val;
Real * stress_val = stress_vect->values;
Real * sigma_dphi_dx_val = sigma_dphi_dx->values;
UInt * elem_filter_val = elem_filter->values;
UInt offset_shapesd_val = spatial_dimension * nb_nodes_per_element;
UInt offset_stress_val = spatial_dimension * spatial_dimension;
UInt offset_sigma_dphi_dx_val = spatial_dimension * nb_nodes_per_element;
for (UInt el = 0; el < nb_element; ++el) {
shapesd_val = shapesd + elem_filter_val[el]*size_of_shapesd*nb_quadrature_points;
for (UInt q = 0; q < nb_quadrature_points; ++q) {
Math::matrix_matrixt(nb_nodes_per_element, spatial_dimension, spatial_dimension,
shapesd_val, stress_val, sigma_dphi_dx_val);
shapesd_val += offset_shapesd_val;
stress_val += offset_stress_val;
sigma_dphi_dx_val += offset_sigma_dphi_dx_val;
}
}
/**
* compute @f$\int \sigma * \frac{\partial \varphi}{\partial X}dX@f$ by @f$ \sum_q \mathbf{B}^t
* \mathbf{\sigma}_q \overline w_q J_q@f$
*/
Vector<Real> * int_sigma_dphi_dx = new Vector<Real>(0,
nb_nodes_per_element * spatial_dimension,
"int_sigma_x_dphi_/_dX");
model->getFEM().integrate(*sigma_dphi_dx, *int_sigma_dphi_dx,
size_of_shapes_derivatives,
*it, ghost_type,
elem_filter);
delete sigma_dphi_dx;
/// assemble
model->getFEM().assembleVector(*int_sigma_dphi_dx, residual,
residual.getNbComponent(),
*it, ghost_type, elem_filter, -1);
delete int_sigma_dphi_dx;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void Material::computePotentialEnergyByElement() {
AKANTU_DEBUG_IN();
const Mesh::ConnectivityTypeList & type_list = model->getFEM().getMesh().getConnectivityTypeList(_not_ghost);
Mesh::ConnectivityTypeList::const_iterator it;
for(it = type_list.begin(); it != type_list.end(); ++it) {
if(model->getFEM().getMesh().getSpatialDimension(*it) != spatial_dimension) continue;
if(potential_energy[*it] == NULL) {
UInt nb_element = element_filter[*it]->getSize();
UInt nb_quadrature_points = FEM::getNbQuadraturePoints(*it);
-
+
std::stringstream sstr; sstr << id << ":potential_energy:"<< *it;
potential_energy[*it] = &(alloc<Real> (sstr.str(), nb_element * nb_quadrature_points,
1,
REAL_INIT_VALUE));
}
computePotentialEnergy(*it);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
Real Material::getPotentialEnergy() {
AKANTU_DEBUG_IN();
Real epot = 0.;
computePotentialEnergyByElement();
/// integrate the potential energy for each type of elements
const Mesh::ConnectivityTypeList & type_list = model->getFEM().getMesh().getConnectivityTypeList(_not_ghost);
Mesh::ConnectivityTypeList::const_iterator it;
for(it = type_list.begin(); it != type_list.end(); ++it) {
if(model->getFEM().getMesh().getSpatialDimension(*it) != spatial_dimension) continue;
epot += model->getFEM().integrate(*potential_energy[*it], *it,
_not_ghost, element_filter[*it]);
}
AKANTU_DEBUG_OUT();
return epot;
}
/* -------------------------------------------------------------------------- */
__END_AKANTU__
diff --git a/model/solid_mechanics/material_parser.hh b/model/solid_mechanics/material_parser.hh
index d4fb9707c..cebc907ff 100644
--- a/model/solid_mechanics/material_parser.hh
+++ b/model/solid_mechanics/material_parser.hh
@@ -1,79 +1,71 @@
/**
* @file material_parser.hh
* @author Guillaume ANCIAUX <anciaux@lsmscluster1.epfl.ch>
* @date Thu Nov 25 11:43:48 2010
*
- * @brief
+ * @brief
*
* @section LICENSE
*
* <insert license here>
*
*/
/* -------------------------------------------------------------------------- */
-
-/* -------------------------------------------------------------------------- */
-
#ifndef __AKANTU_MATERIAL_PARSER_HH__
#define __AKANTU_MATERIAL_PARSER_HH__
__BEGIN_AKANTU__
class MaterialParser {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
-
- MaterialParser(){};
- virtual ~MaterialParser(){infile.close();};
-
+
+ MaterialParser() : current_line(0) {};
+ virtual ~MaterialParser(){ infile.close(); };
+
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// open a file to parse
void open(const std::string & filename);
/// read the file and return the next material type
std::string getNextMaterialType();
/// read properties and instanciate a given material object
- template <typename M>
- Material * readMaterialObject(SolidMechanicsModel & model,MaterialID & mat_id);
+ template <typename Mat>
+ Material * readMaterialObject(SolidMechanicsModel & model, MaterialID & mat_id);
-
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
-
+
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
+ inline void my_getline();
-
- inline void my_getline() {
- std::getline(infile, line); //read the line
- if (!(infile.flags() & (std::ios::failbit | std::ios::eofbit)))
- ++current_line;
- size_t pos = line.find("#"); //remove the comment
- line = line.substr(0, pos);
- trim(line); // remove unnecessary spaces
- }
-
+ /// number of the current line
UInt current_line;
+
+ /// material file
std::ifstream infile;
+
+ /// current read line
std::string line;
};
#include "material_parser_inline_impl.cc"
__END_AKANTU__
#endif
diff --git a/model/solid_mechanics/material_parser_inline_impl.cc b/model/solid_mechanics/material_parser_inline_impl.cc
index 83fe33b72..5eed96361 100644
--- a/model/solid_mechanics/material_parser_inline_impl.cc
+++ b/model/solid_mechanics/material_parser_inline_impl.cc
@@ -1,83 +1,93 @@
/**
* @file material_parser_inline_impl.cc
- * @author Guillaume ANCIAUX <anciaux@lsmscluster1.epfl.ch>
+ * @author Guillaume ANCIAUX <guillaume.anciaux@epfl.ch>
+ * @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Fri Nov 26 08:32:38 2010
*
- * @brief
+ * @brief
*
* @section LICENSE
*
* <insert license here>
*
*/
/* -------------------------------------------------------------------------- */
-
-template <typename M>
+template <typename Mat>
inline Material * MaterialParser::readMaterialObject(SolidMechanicsModel & model,MaterialID & mat_id){
std::string keyword;
std::string value;
-
+
/// instanciate the material object
- Material * mat = new M(model, mat_id);
+ Material * mat = new Mat(model, mat_id);
/// read the material properties
my_getline();
-
+
while(line[0] != ']') {
size_t pos = line.find("=");
if(pos == std::string::npos)
AKANTU_DEBUG_ERROR("Malformed material file : line must be \"key = value\" at line"
<< current_line);
-
+
keyword = line.substr(0, pos); trim(keyword);
value = line.substr(pos + 1); trim(value);
-
+
try {
mat->setParam(keyword, value, mat_id);
- } catch (Exception ex) {
+ } catch (debug::Exception ex) {
AKANTU_DEBUG_ERROR("Malformed material file : error in setParam \""
<< ex.info() << "\" at line " << current_line);
}
-
+
my_getline();
}
return mat;
}
+/* -------------------------------------------------------------------------- */
inline std::string MaterialParser::getNextMaterialType(){
while(infile.good()) {
my_getline();
-
+
// if empty line continue
if(line.empty()) continue;
-
+
std::stringstream sstr(line);
std::string keyword;
std::string value;
-
+
sstr >> keyword;
to_lower(keyword);
/// if found a material deccription then stop
- /// and prepare the things for further reading
+ /// and prepare the things for further reading
if(keyword == "material") {
std::string type; sstr >> type;
to_lower(type);
std::string obracket; sstr >> obracket;
if(obracket != "[")
AKANTU_DEBUG_ERROR("Malformed material file : missing [ at line " << current_line);
return type;
- }
+ }
}
return "";
}
+/* -------------------------------------------------------------------------- */
+inline void MaterialParser::my_getline() {
+ std::getline(infile, line); //read the line
+ if (!(infile.flags() & (std::ios::failbit | std::ios::eofbit)))
+ ++current_line;
+ size_t pos = line.find("#"); //remove the comment
+ line = line.substr(0, pos);
+ trim(line); // remove unnecessary spaces
+}
+/* -------------------------------------------------------------------------- */
inline void MaterialParser::open(const std::string & filename){
infile.open(filename.c_str());
current_line = 0;
-
+
if(!infile.good()) {
AKANTU_DEBUG_ERROR("Cannot open file " << filename);
}
}
-
diff --git a/model/solid_mechanics/materials/material_elastic_inline_impl.cc b/model/solid_mechanics/materials/material_elastic_inline_impl.cc
index e2672eb83..381152500 100644
--- a/model/solid_mechanics/materials/material_elastic_inline_impl.cc
+++ b/model/solid_mechanics/materials/material_elastic_inline_impl.cc
@@ -1,48 +1,49 @@
/**
* @file material_elastic_inline_impl.cc
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Tue Jul 27 11:57:43 2010
*
* @brief Implementation of the inline functions of the material elastic
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
inline void MaterialElastic::computeStress(Real * F, Real * sigma) {
Real trace = F[0] + F[4] + F[8]; /// \F_{11} + \F_{22} + \F_{33}
/// \sigma_{ij} = \lamda * \F_{kk} * \delta_{ij} + 2 * \mu * \F_{ij}
sigma[0] = lambda * trace + 2*mu*F[0];
sigma[4] = lambda * trace + 2*mu*F[4];
sigma[8] = lambda * trace + 2*mu*F[8];
sigma[1] = sigma[3] = mu * (F[1] + F[3]);
sigma[2] = sigma[6] = mu * (F[2] + F[6]);
sigma[5] = sigma[7] = mu * (F[5] + F[7]);
}
/* -------------------------------------------------------------------------- */
inline void MaterialElastic::computePotentialEnergy(Real * F, Real * sigma, Real * epot) {
*epot = 0.;
for (UInt i = 0, t = 0; i < spatial_dimension; ++i)
for (UInt j = 0; j < spatial_dimension; ++j, ++t)
- *epot += sigma[t] * (F[t] - (i == j));
+ (*epot) += sigma[t] * (F[t] - (i == j));
+
*epot *= .5;
}
/* -------------------------------------------------------------------------- */
inline Real MaterialElastic::celerity() {
return sqrt(E/rho);
}
/* -------------------------------------------------------------------------- */
inline Real MaterialElastic::getStableTimeStep(Real h) {
return (h/celerity());
}
diff --git a/model/solid_mechanics/solid_mechanics_model.cc b/model/solid_mechanics/solid_mechanics_model.cc
index 5560182ae..297e5f166 100644
--- a/model/solid_mechanics/solid_mechanics_model.cc
+++ b/model/solid_mechanics/solid_mechanics_model.cc
@@ -1,433 +1,455 @@
/**
* @file solid_mechanics_model.cc
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Thu Jul 22 14:35:38 2010
*
* @brief Implementation of the SolidMechanicsModel class
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -------------------------------------------------------------------------- */
#include "solid_mechanics_model.hh"
#include "aka_math.hh"
#include "integration_scheme_2nd_order.hh"
#include "static_communicator.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh,
UInt spatial_dimension,
const ModelID & id,
const MemoryID & memory_id) :
Model(mesh, spatial_dimension, id, memory_id),
time_step(NAN), f_m2a(1.0),
- integrator(new CentralDifference()) {
+ integrator(new CentralDifference()),
+ increment_flag(false) {
AKANTU_DEBUG_IN();
this->displacement = NULL;
this->mass = NULL;
this->velocity = NULL;
this->acceleration = NULL;
this->force = NULL;
this->residual = NULL;
this->boundary = NULL;
this->increment = NULL;
for(UInt t = _not_defined; t < _max_element_type; ++t) {
this->element_material[t] = NULL;
this->ghost_element_material[t] = NULL;
}
registerTag(_gst_smm_mass, "Mass");
registerTag(_gst_smm_residual, "Explicit Residual");
registerTag(_gst_smm_boundary, "Boundary conditions");
materials.clear();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
SolidMechanicsModel::~SolidMechanicsModel() {
AKANTU_DEBUG_IN();
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
delete *mat_it;
}
materials.clear();
delete integrator;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initVectors() {
AKANTU_DEBUG_IN();
UInt nb_nodes = fem->getMesh().getNbNodes();
std::stringstream sstr_disp; sstr_disp << id << ":displacement";
std::stringstream sstr_mass; sstr_mass << id << ":mass";
std::stringstream sstr_velo; sstr_velo << id << ":velocity";
std::stringstream sstr_acce; sstr_acce << id << ":acceleration";
std::stringstream sstr_forc; sstr_forc << id << ":force";
std::stringstream sstr_resi; sstr_resi << id << ":residual";
std::stringstream sstr_boun; sstr_boun << id << ":boundary";
displacement = &(alloc<Real>(sstr_disp.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE));
mass = &(alloc<Real>(sstr_mass.str(), nb_nodes, 1)); // \todo see if it must not be spatial_dimension
velocity = &(alloc<Real>(sstr_velo.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE));
acceleration = &(alloc<Real>(sstr_acce.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE));
force = &(alloc<Real>(sstr_forc.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE));
residual = &(alloc<Real>(sstr_resi.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE));
boundary = &(alloc<bool>(sstr_boun.str(), nb_nodes, spatial_dimension, false));
std::stringstream sstr_curp; sstr_curp << id << ":current_position_tmp";
current_position = &(alloc<Real>(sstr_curp.str(), 0, spatial_dimension, REAL_INIT_VALUE));
const Mesh::ConnectivityTypeList & type_list = fem->getMesh().getConnectivityTypeList();
Mesh::ConnectivityTypeList::const_iterator it;
for(it = type_list.begin(); it != type_list.end(); ++it) {
UInt nb_element = fem->getMesh().getNbElement(*it);
if(!element_material[*it]) {
std::stringstream sstr_elma; sstr_elma << id << ":element_material:" << *it;
element_material[*it] = &(alloc<UInt>(sstr_elma.str(), nb_element, 1, 0));
}
}
const Mesh::ConnectivityTypeList & ghost_type_list =
fem->getMesh().getConnectivityTypeList(_ghost);
for(it = ghost_type_list.begin(); it != ghost_type_list.end(); ++it) {
UInt nb_element = fem->getMesh().getNbGhostElement(*it);
std::stringstream sstr_elma; sstr_elma << id << ":ghost_element_material:" << *it;
ghost_element_material[*it] = &(alloc<UInt>(sstr_elma.str(), nb_element, 1, 0));
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initModel() {
fem->initShapeFunctions(_not_ghost);
fem->initShapeFunctions(_ghost);
}
/* -------------------------------------------------------------------------- */
-void SolidMechanicsModel::updateResidual() {
- AKANTU_DEBUG_IN();
-
+void SolidMechanicsModel::updateCurrentPosition() {
UInt nb_nodes = fem->getMesh().getNbNodes();
- residual->resize(nb_nodes);
current_position->resize(nb_nodes);
//Vector<Real> * current_position = new Vector<Real>(nb_nodes, spatial_dimension, NAN, "position");
Real * current_position_val = current_position->values;
Real * position_val = fem->getMesh().getNodes().values;
Real * displacement_val = displacement->values;
/// compute current_position = initial_position + displacement
memcpy(current_position_val, position_val, nb_nodes*spatial_dimension*sizeof(Real));
for (UInt n = 0; n < nb_nodes*spatial_dimension; ++n) {
*current_position_val++ += *displacement_val++;
}
+}
+
+/* -------------------------------------------------------------------------- */
+void SolidMechanicsModel::updateResidual(bool need_update_current_position) {
+ AKANTU_DEBUG_IN();
+
+ UInt nb_nodes = fem->getMesh().getNbNodes();
+
+ residual->resize(nb_nodes);
+
+ if (need_update_current_position) updateCurrentPosition();
/// start synchronization
asynchronousSynchronize(_gst_smm_residual);
/// copy the forces in residual for boundary conditions
memcpy(residual->values, force->values, nb_nodes*spatial_dimension*sizeof(Real));
/// call update residual on each local elements
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
(*mat_it)->updateResidual(*current_position, _not_ghost);
}
/// finalize communications
waitEndSynchronize(_gst_smm_residual);
/// call update residual on each ghost elements
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
(*mat_it)->updateResidual(*current_position, _ghost);
}
// current_position;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateAcceleration() {
AKANTU_DEBUG_IN();
UInt nb_nodes = acceleration->getSize();
UInt nb_degre_of_freedom = acceleration->getNbComponent();
Real * mass_val = mass->values;
Real * residual_val = residual->values;
Real * accel_val = acceleration->values;
bool * boundary_val = boundary->values;
for (UInt n = 0; n < nb_nodes; ++n) {
for (UInt d = 0; d < nb_degre_of_freedom; d++) {
if(!(*boundary_val)) {
*accel_val = f_m2a * *residual_val / *mass_val;
}
residual_val++;
accel_val++;
boundary_val++;
}
mass_val++;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::explicitPred() {
AKANTU_DEBUG_IN();
if(increment_flag) {
memcpy(increment->values, displacement->values, displacement->getSize()*displacement->getNbComponent()*sizeof(Real));
}
integrator->integrationSchemePred(time_step,
*displacement,
*velocity,
*acceleration,
*boundary);
if(increment_flag) {
Real * inc_val = increment->values;
Real * dis_val = displacement->values;
UInt nb_nodes = displacement->getSize();
for (UInt n = 0; n < nb_nodes; ++n) {
*inc_val = *dis_val - *inc_val;
*inc_val++;
*dis_val++;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::explicitCorr() {
AKANTU_DEBUG_IN();
integrator->integrationSchemeCorr(time_step,
*displacement,
*velocity,
*acceleration,
*boundary);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::synchronizeBoundaries() {
AKANTU_DEBUG_IN();
synchronize(_gst_smm_boundary);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::setIncrementFlagOn() {
AKANTU_DEBUG_IN();
if(!increment) {
UInt nb_nodes = fem->getMesh().getNbNodes();
std::stringstream sstr_inc; sstr_inc << id << ":increment";
increment = &(alloc<Real>(sstr_inc.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE));
}
increment_flag = true;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getStableTimeStep() {
AKANTU_DEBUG_IN();
Material ** mat_val = &(materials.at(0));
Real min_dt = std::numeric_limits<Real>::max();
Real * coord = fem->getMesh().getNodes().values;
Real * disp_val = displacement->values;
const Mesh::ConnectivityTypeList & type_list = fem->getMesh().getConnectivityTypeList();
Mesh::ConnectivityTypeList::const_iterator it;
for(it = type_list.begin(); it != type_list.end(); ++it) {
if(fem->getMesh().getSpatialDimension(*it) != spatial_dimension) continue;
UInt nb_nodes_per_element = fem->getMesh().getNbNodesPerElement(*it);
UInt nb_element = fem->getMesh().getNbElement(*it);
UInt * conn = fem->getMesh().getConnectivity(*it).values;
UInt * elem_mat_val = element_material[*it]->values;
Real * u = new Real[nb_nodes_per_element*spatial_dimension];
for (UInt el = 0; el < nb_element; ++el) {
UInt el_offset = el * nb_nodes_per_element;
for (UInt n = 0; n < nb_nodes_per_element; ++n) {
UInt offset_conn = conn[el_offset + n] * spatial_dimension;
memcpy(u + n * spatial_dimension,
coord + offset_conn,
spatial_dimension * sizeof(Real));
for (UInt i = 0; i < spatial_dimension; ++i) {
u[n * spatial_dimension + i] += disp_val[offset_conn + i];
}
}
Real el_size = fem->getElementInradius(u, *it);
Real el_dt = mat_val[elem_mat_val[el]]->getStableTimeStep(el_size);
min_dt = min_dt > el_dt ? el_dt : min_dt;
}
delete [] u;
}
/// reduction min over all processors
allReduce(&min_dt, _so_min);
AKANTU_DEBUG_OUT();
return min_dt;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getPotentialEnergy() {
AKANTU_DEBUG_IN();
Real epot = 0.;
/// call update residual on each local elements
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
epot += (*mat_it)->getPotentialEnergy();
}
/// reduction sum over all processors
allReduce(&epot, _so_sum);
AKANTU_DEBUG_OUT();
return epot;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getKineticEnergy() {
AKANTU_DEBUG_IN();
Real ekin = 0.;
UInt nb_nodes = fem->getMesh().getNbNodes();
- Vector<Real> * v_square = new Vector<Real>(nb_nodes, 1, "v_square");
+ // Vector<Real> * v_square = new Vector<Real>(nb_nodes, 1, "v_square");
Real * vel_val = velocity->values;
- Real * v_s_val = v_square->values;
+ Real * mass_val = mass->values;
for (UInt n = 0; n < nb_nodes; ++n) {
- *v_s_val = 0;
- for (UInt s = 0; s < spatial_dimension; ++s) {
- *v_s_val += *vel_val * *vel_val;
+ Real v2 = 0;
+ for (UInt i = 0; i < spatial_dimension; ++i) {
+ v2 += *vel_val * *vel_val;
vel_val++;
}
- v_s_val++;
+ ekin += *mass_val * v2;
+ mass_val++;
}
- Material ** mat_val = &(materials.at(0));
- const Mesh:: ConnectivityTypeList & type_list = fem->getMesh().getConnectivityTypeList();
- Mesh::ConnectivityTypeList::const_iterator it;
- for(it = type_list.begin(); it != type_list.end(); ++it) {
- if(fem->getMesh().getSpatialDimension(*it) != spatial_dimension) continue;
+ // Real * v_s_val = v_square->values;
- UInt nb_quadrature_points = FEM::getNbQuadraturePoints(*it);
- UInt nb_element = fem->getMesh().getNbElement(*it);
+ // for (UInt n = 0; n < nb_nodes; ++n) {
+ // *v_s_val = 0;
+ // for (UInt s = 0; s < spatial_dimension; ++s) {
+ // *v_s_val += *vel_val * *vel_val;
+ // vel_val++;
+ // }
+ // v_s_val++;
+ // }
- Vector<Real> * v_square_el = new Vector<Real>(nb_element*nb_quadrature_points, 1, "v_square per element");
+ // Material ** mat_val = &(materials.at(0));
- fem->interpolateOnQuadraturePoints(*v_square, *v_square_el, 1, *it);
+ // const Mesh:: ConnectivityTypeList & type_list = fem->getMesh().getConnectivityTypeList();
+ // Mesh::ConnectivityTypeList::const_iterator it;
+ // for(it = type_list.begin(); it != type_list.end(); ++it) {
+ // if(fem->getMesh().getSpatialDimension(*it) != spatial_dimension) continue;
- Real * v_square_el_val = v_square_el->values;
- UInt * elem_mat_val = element_material[*it]->values;
+ // UInt nb_quadrature_points = FEM::getNbQuadraturePoints(*it);
+ // UInt nb_element = fem->getMesh().getNbElement(*it);
- for (UInt el = 0; el < nb_element; ++el) {
- Real rho = mat_val[elem_mat_val[el]]->getRho();
- for (UInt q = 0; q < nb_quadrature_points; ++q) {
- *v_square_el_val *= rho;
- v_square_el_val++;
- }
- }
+ // Vector<Real> * v_square_el = new Vector<Real>(nb_element * nb_quadrature_points, 1, "v_square per element");
- ekin += fem->integrate(*v_square_el, *it);
+ // fem->interpolateOnQuadraturePoints(*v_square, *v_square_el, 1, *it);
- delete v_square_el;
- }
- delete v_square;
+ // Real * v_square_el_val = v_square_el->values;
+ // UInt * elem_mat_val = element_material[*it]->values;
+
+ // for (UInt el = 0; el < nb_element; ++el) {
+ // Real rho = mat_val[elem_mat_val[el]]->getRho();
+ // for (UInt q = 0; q < nb_quadrature_points; ++q) {
+ // *v_square_el_val *= rho;
+ // v_square_el_val++;
+ // }
+ // }
+
+ // ekin += fem->integrate(*v_square_el, *it);
+
+ // delete v_square_el;
+ // }
+ // delete v_square;
/// reduction sum over all processors
allReduce(&ekin, _so_sum);
AKANTU_DEBUG_OUT();
return ekin * .5;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::printself(std::ostream & stream, int indent) const {
std::string space;
for(Int i = 0; i < indent; i++, space += AKANTU_INDENT);
stream << space << "Solid Mechanics Model [" << std::endl;
stream << space << " + id : " << id << std::endl;
stream << space << " + spatial dimension : " << spatial_dimension << std::endl;
stream << space << " + fem [" << std::endl;
fem->printself(stream, indent + 2);
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << " + nodals information [" << std::endl;
displacement->printself(stream, indent + 2);
mass ->printself(stream, indent + 2);
velocity ->printself(stream, indent + 2);
acceleration->printself(stream, indent + 2);
force ->printself(stream, indent + 2);
residual ->printself(stream, indent + 2);
boundary ->printself(stream, indent + 2);
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << " + connectivity type information [" << std::endl;
const Mesh::ConnectivityTypeList & type_list = fem->getMesh().getConnectivityTypeList();
Mesh::ConnectivityTypeList::const_iterator it;
for(it = type_list.begin(); it != type_list.end(); ++it) {
stream << space << AKANTU_INDENT << AKANTU_INDENT << " + " << *it <<" [" << std::endl;
element_material[*it]->printself(stream, indent + 3);
stream << space << AKANTU_INDENT << AKANTU_INDENT << "]" << std::endl;
}
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
__END_AKANTU__
diff --git a/model/solid_mechanics/solid_mechanics_model.hh b/model/solid_mechanics/solid_mechanics_model.hh
index ea28c4efd..df8a1e11e 100644
--- a/model/solid_mechanics/solid_mechanics_model.hh
+++ b/model/solid_mechanics/solid_mechanics_model.hh
@@ -1,256 +1,259 @@
/**
* @file solid_mechanics_model.hh
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date[creation] Thu Jul 22 11:51:06 2010
* @date[last modification] Thu Oct 14 14:00:06 2010
*
* @brief Model of Solid Mechanics
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_SOLID_MECHANICS_MODEL_HH__
#define __AKANTU_SOLID_MECHANICS_MODEL_HH__
/* -------------------------------------------------------------------------- */
#include <fstream>
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "model.hh"
#include "material.hh"
#include "material_parser.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
// class Material;
class IntegrationScheme2ndOrder;
class Contact;
}
__BEGIN_AKANTU__
class SolidMechanicsModel : public Model {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
SolidMechanicsModel(UInt spatial_dimension,
const ModelID & id = "solid_mechanics_model",
const MemoryID & memory_id = 0);
SolidMechanicsModel(Mesh & mesh,
UInt spatial_dimension = 0,
const ModelID & id = "solid_mechanics_model",
const MemoryID & memory_id = 0);
virtual ~SolidMechanicsModel();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// allocate all vectors
void initVectors();
/// read the material files to instantiate all the materials
void readMaterials(const std::string & filename);
/// read a custom material with a keyword and class as template
template <typename M>
UInt readCustomMaterial(const std::string & filename,
const std::string & keyword);
/// read properties part of a material file and create the material
template <typename M>
Material * readMaterialProperties(std::ifstream & infile,
MaterialID mat_id,
UInt ¤t_line);
/// initialize all internal arrays for materials
void initMaterials();
/// initialize the model
void initModel();
/// assemble the lumped mass matrix
void assembleMassLumped();
+ /// update the current position vector
+ void updateCurrentPosition();
+
/// assemble the residual for the explicit scheme
- void updateResidual();
+ void updateResidual(bool need_update_current_position = true);
/// compute the acceleration from the residual
void updateAcceleration();
/// explicit integration predictor
void explicitPred();
/// explicit integration corrector
void explicitCorr();
/// synchronize the ghost element boundaries values
void synchronizeBoundaries();
/// function to print the containt of the class
virtual void printself(std::ostream & stream, int indent = 0) const;
/// integrate a force on the boundary by providing a stress tensor
void computeForcesByStressTensor(const Vector<Real> & stresses, const ElementType & type);
/// integrate a force on the boundary by providing a traction vector
void computeForcesByTractionVector(const Vector<Real> & tractions, const ElementType & type);
/// compute force vector from a function(x,y,z) that describe stresses
void computeForcesFromFunction(void (*myf)(double *,double *), BoundaryFunctionType function_type);
private:
/// assemble the lumped mass matrix for local and ghost elements
void assembleMassLumped(GhostType ghost_type);
/// assemble the lumped mass matrix for local and ghost elements
void assembleMassLumpedRowSum(GhostType ghost_type, const ElementType type);
/// assemble the lumped mass matrix for local and ghost elements
void assembleMassLumpedDiagonalScaling(GhostType ghost_type, const ElementType type);
/* ------------------------------------------------------------------------ */
/* Ghost Synchronizer inherited members */
/* ------------------------------------------------------------------------ */
public:
inline virtual UInt getNbDataToPack(const Element & element,
GhostSynchronizationTag tag) const;
inline virtual UInt getNbDataToUnpack(const Element & element,
GhostSynchronizationTag tag) const;
inline virtual void packData(Real ** buffer,
const Element & element,
GhostSynchronizationTag tag) const;
inline virtual void unpackData(Real ** buffer,
const Element & element,
GhostSynchronizationTag tag) const;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
AKANTU_GET_MACRO(TimeStep, time_step, Real);
AKANTU_SET_MACRO(TimeStep, time_step, Real);
AKANTU_GET_MACRO(F_M2A, f_m2a, Real);
AKANTU_SET_MACRO(F_M2A, f_m2a, Real);
AKANTU_GET_MACRO(Displacement, *displacement, Vector<Real> &);
AKANTU_GET_MACRO(CurrentPosition, *current_position, const Vector<Real> &);
AKANTU_GET_MACRO(Increment, *increment, const Vector<Real> &);
AKANTU_GET_MACRO(Mass, *mass, const Vector<Real> &);
AKANTU_GET_MACRO(Velocity, *velocity, Vector<Real> &);
AKANTU_GET_MACRO(Acceleration, *acceleration, Vector<Real> &);
AKANTU_GET_MACRO(Force, *force, Vector<Real> &);
AKANTU_GET_MACRO(Residual, *residual, const Vector<Real> &);
AKANTU_GET_MACRO(Boundary, *boundary, Vector<bool> &);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementMaterial, element_material, Vector<UInt> &);
inline Material & getMaterial(UInt mat_index);
/// compute the stable time step
Real getStableTimeStep();
// void setPotentialEnergyFlagOn();
// void setPotentialEnergyFlagOff();
Real getPotentialEnergy();
Real getKineticEnergy();
AKANTU_SET_MACRO(Contact, contact, Contact *);
void setIncrementFlagOn();
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
/// time step
Real time_step;
/// conversion coefficient form force/mass to acceleration
Real f_m2a;
/// displacements array
Vector<Real> * displacement;
/// lumped mass array
Vector<Real> * mass;
/// velocities array
Vector<Real> * velocity;
/// accelerations array
Vector<Real> * acceleration;
/// forces array
Vector<Real> * force;
/// residuals array
Vector<Real> * residual;
/// boundaries array
Vector<bool> * boundary;
/// array of current position used during update residual
Vector<Real> * current_position;
/// materials of all elements
ByElementTypeUInt element_material;
/// materials of all ghost elements
ByElementTypeUInt ghost_element_material;
/// list of used materials
std::vector<Material *> materials;
/// integration scheme of second order used
IntegrationScheme2ndOrder * integrator;
/// increment of displacement
Vector<Real> * increment;
/// flag defining if the increment must be computed or not
bool increment_flag;
/// object to resolve the contact
Contact * contact;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#include "solid_mechanics_model_inline_impl.cc"
/// standard output stream operator
inline std::ostream & operator <<(std::ostream & stream, const SolidMechanicsModel & _this)
{
_this.printself(stream);
return stream;
}
__END_AKANTU__
#endif /* __AKANTU_SOLID_MECHANICS_MODEL_HH__ */
diff --git a/synchronizer/static_communicator_mpi.hh b/synchronizer/static_communicator_mpi.hh
index 16260eee8..56a8a649f 100644
--- a/synchronizer/static_communicator_mpi.hh
+++ b/synchronizer/static_communicator_mpi.hh
@@ -1,110 +1,111 @@
/**
* @file static_communicator_mpi.hh
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Thu Sep 2 19:59:58 2010
*
* @brief class handling parallel communication trough MPI
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_STATIC_COMMUNICATOR_MPI_HH__
#define __AKANTU_STATIC_COMMUNICATOR_MPI_HH__
/* -------------------------------------------------------------------------- */
#include <mpi.h>
/* -------------------------------------------------------------------------- */
#include "static_communicator.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
class CommunicationRequestMPI : public CommunicationRequest {
public:
inline CommunicationRequestMPI(UInt source, UInt dest);
inline ~CommunicationRequestMPI();
inline MPI_Request * getMPIRequest() { return request; };
private:
MPI_Request * request;
};
/* -------------------------------------------------------------------------- */
class StaticCommunicatorMPI : public virtual StaticCommunicator {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
inline StaticCommunicatorMPI(int * argc, char *** argv);
inline virtual ~StaticCommunicatorMPI();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
inline void send(UInt * buffer, Int size, Int receiver, Int tag);
inline void send(Real * buffer, Int size, Int receiver, Int tag);
inline void receive(UInt * buffer, Int size, Int sender, Int tag);
inline void receive(Real * buffer, Int size, Int sender, Int tag);
inline CommunicationRequest * asyncSend(UInt * buffer, Int size, Int receiver, Int tag);
inline CommunicationRequest * asyncSend(Real * buffer, Int size, Int receiver, Int tag);
inline CommunicationRequest * asyncReceive(UInt * buffer, Int size,
Int sender, Int tag);
inline CommunicationRequest * asyncReceive(Real * buffer, Int size,
Int sender, Int tag);
inline bool testRequest(CommunicationRequest * request);
inline void wait(CommunicationRequest * request);
inline void waitAll(std::vector<CommunicationRequest *> & requests);
inline void barrier();
inline void allReduce(Real * values, UInt nb_values, const SynchronizerOperation & op);
inline void allReduce(UInt * values, UInt nb_values, const SynchronizerOperation & op);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
- inline void setMPIComm(MPI_Comm comm);
+ inline void setMPICommunicator(MPI_Comm comm);
+ inline MPI_Comm getMPICommunicator();
inline Int getNbProc() { return psize; };
inline Int whoAmI() { return prank; };
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
MPI_Comm communicator;
static MPI_Op synchronizer_operation_to_mpi_op[_so_null + 1];
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#include "static_communicator_mpi_inline_impl.cc"
__END_AKANTU__
#endif /* __AKANTU_STATIC_COMMUNICATOR_MPI_HH__ */
diff --git a/synchronizer/static_communicator_mpi_inline_impl.cc b/synchronizer/static_communicator_mpi_inline_impl.cc
index 81002f33a..bfbf6f94f 100644
--- a/synchronizer/static_communicator_mpi_inline_impl.cc
+++ b/synchronizer/static_communicator_mpi_inline_impl.cc
@@ -1,163 +1,226 @@
/**
- * @file static_communicator_mpi_inline_impl.hh
+ * @file static_communicator_mpi_inline_impl.cc
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Thu Sep 2 15:10:51 2010
*
* @brief implementation of the mpi communicator
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -------------------------------------------------------------------------- */
inline CommunicationRequestMPI::CommunicationRequestMPI(UInt source, UInt dest) :
CommunicationRequest(source, dest) {
request = new MPI_Request;
}
/* -------------------------------------------------------------------------- */
inline CommunicationRequestMPI::~CommunicationRequestMPI() {
delete request;
}
/* -------------------------------------------------------------------------- */
inline StaticCommunicatorMPI::StaticCommunicatorMPI(int * argc, char *** argv) {
MPI_Init(argc, argv);
- setMPIComm(MPI_COMM_WORLD);
+ setMPICommunicator(MPI_COMM_WORLD);
}
/* -------------------------------------------------------------------------- */
inline StaticCommunicatorMPI::~StaticCommunicatorMPI() {
MPI_Finalize();
}
/* -------------------------------------------------------------------------- */
-inline void StaticCommunicatorMPI::setMPIComm(MPI_Comm comm) {
+inline void StaticCommunicatorMPI::setMPICommunicator(MPI_Comm comm) {
communicator = comm;
MPI_Comm_rank(communicator, &prank);
MPI_Comm_size(communicator, &psize);
}
+/* -------------------------------------------------------------------------- */
+inline MPI_Comm StaticCommunicatorMPI::getMPICommunicator() {
+ return communicator;
+}
+
/* -------------------------------------------------------------------------- */
inline void StaticCommunicatorMPI::send(UInt * buffer, Int size,
Int receiver, Int tag) {
- int ret = MPI_Send(buffer, size, MPI_UNSIGNED, receiver, tag, communicator);
+#if !defined(AKANTU_NDEBUG)
+ int ret =
+#endif
+ MPI_Send(buffer, size, MPI_UNSIGNED, receiver, tag, communicator);
+
AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Send.");
}
/* -------------------------------------------------------------------------- */
inline void StaticCommunicatorMPI::send(Real * buffer, Int size,
Int receiver, Int tag) {
- int ret = MPI_Send(buffer, size, MPI_DOUBLE, receiver, tag, communicator);
+#if !defined(AKANTU_NDEBUG)
+ int ret =
+#endif
+ MPI_Send(buffer, size, MPI_DOUBLE, receiver, tag, communicator);
+
AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Send.");
}
/* -------------------------------------------------------------------------- */
inline void StaticCommunicatorMPI::receive(UInt * buffer, Int size,
Int sender, Int tag) {
MPI_Status status;
- int ret = MPI_Recv(buffer, size, MPI_UNSIGNED, sender, tag, communicator, &status);
+
+#if !defined(AKANTU_NDEBUG)
+ int ret =
+#endif
+ MPI_Recv(buffer, size, MPI_UNSIGNED, sender, tag, communicator, &status);
+
AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Recv.");
}
/* -------------------------------------------------------------------------- */
inline void StaticCommunicatorMPI::receive(Real * buffer, Int size,
Int sender, Int tag) {
MPI_Status status;
- int ret = MPI_Recv(buffer, size, MPI_DOUBLE, sender, tag, communicator, &status);
+
+#if !defined(AKANTU_NDEBUG)
+ int ret =
+#endif
+ MPI_Recv(buffer, size, MPI_DOUBLE, sender, tag, communicator, &status);
+
AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Recv.");
}
/* -------------------------------------------------------------------------- */
inline CommunicationRequest * StaticCommunicatorMPI::asyncSend(UInt * buffer, Int size,
Int receiver, Int tag) {
CommunicationRequestMPI * request = new CommunicationRequestMPI(prank, receiver);
- int ret = MPI_Isend(buffer, size, MPI_UNSIGNED, receiver, tag, communicator, request->getMPIRequest());
+
+#if !defined(AKANTU_NDEBUG)
+ int ret =
+#endif
+ MPI_Isend(buffer, size, MPI_UNSIGNED, receiver, tag, communicator, request->getMPIRequest());
+
AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Isend.");
return request;
}
/* -------------------------------------------------------------------------- */
inline CommunicationRequest * StaticCommunicatorMPI::asyncSend(Real * buffer, Int size,
Int receiver, Int tag) {
CommunicationRequestMPI * request = new CommunicationRequestMPI(prank, receiver);
- int ret = MPI_Isend(buffer, size, MPI_DOUBLE, receiver, tag, communicator, request->getMPIRequest());
+
+#if !defined(AKANTU_NDEBUG)
+ int ret =
+#endif
+ MPI_Isend(buffer, size, MPI_DOUBLE, receiver, tag, communicator, request->getMPIRequest());
+
AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Isend.");
return request;
}
/* -------------------------------------------------------------------------- */
inline CommunicationRequest * StaticCommunicatorMPI::asyncReceive(UInt * buffer, Int size,
Int sender, Int tag) {
CommunicationRequestMPI * request = new CommunicationRequestMPI(sender, prank);
- int ret = MPI_Irecv(buffer, size, MPI_UNSIGNED, sender, tag, communicator, request->getMPIRequest());
+
+#if !defined(AKANTU_NDEBUG)
+ int ret =
+#endif
+ MPI_Irecv(buffer, size, MPI_UNSIGNED, sender, tag, communicator, request->getMPIRequest());
+
AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Irecv.");
return request;
}
/* -------------------------------------------------------------------------- */
inline CommunicationRequest * StaticCommunicatorMPI::asyncReceive(Real * buffer, Int size,
Int sender, Int tag) {
CommunicationRequestMPI * request = new CommunicationRequestMPI(sender, prank);
- int ret = MPI_Irecv(buffer, size, MPI_DOUBLE, sender, tag, communicator, request->getMPIRequest());
+
+#if !defined(AKANTU_NDEBUG)
+ int ret =
+#endif
+ MPI_Irecv(buffer, size, MPI_DOUBLE, sender, tag, communicator, request->getMPIRequest());
+
AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Irecv.");
return request;
}
/* -------------------------------------------------------------------------- */
inline bool StaticCommunicatorMPI::testRequest(CommunicationRequest * request) {
MPI_Status status;
int flag;
CommunicationRequestMPI * req_mpi = static_cast<CommunicationRequestMPI *>(request);
MPI_Request * req = req_mpi->getMPIRequest();
- int ret = MPI_Test(req, &flag, &status);
+#if !defined(AKANTU_NDEBUG)
+ int ret =
+#endif
+
+ MPI_Test(req, &flag, &status);
+
AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Test.");
return (flag != 0);
}
/* -------------------------------------------------------------------------- */
inline void StaticCommunicatorMPI::wait(CommunicationRequest * request) {
MPI_Status status;
CommunicationRequestMPI * req_mpi = static_cast<CommunicationRequestMPI *>(request);
MPI_Request * req = req_mpi->getMPIRequest();
- int ret = MPI_Wait(req, &status);
+#if !defined(AKANTU_NDEBUG)
+ int ret =
+#endif
+ MPI_Wait(req, &status);
+
AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Wait.");
}
/* -------------------------------------------------------------------------- */
inline void StaticCommunicatorMPI::waitAll(std::vector<CommunicationRequest *> & requests) {
MPI_Status status;
std::vector<CommunicationRequest *>::iterator it;
for(it = requests.begin(); it != requests.end(); ++it) {
MPI_Request * req = static_cast<CommunicationRequestMPI *>(*it)->getMPIRequest();
- int ret = MPI_Wait(req, &status);
+
+#if !defined(AKANTU_NDEBUG)
+ int ret =
+#endif
+ MPI_Wait(req, &status);
+
AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Wait.");
}
}
/* -------------------------------------------------------------------------- */
inline void StaticCommunicatorMPI::barrier() {
MPI_Barrier(communicator);
}
/* -------------------------------------------------------------------------- */
inline void StaticCommunicatorMPI::allReduce(UInt * values, UInt nb_values,
const SynchronizerOperation & op) {
- int ret = MPI_Allreduce(MPI_IN_PLACE, values, nb_values, MPI_UNSIGNED,
- synchronizer_operation_to_mpi_op[op],
- communicator);
+#if !defined(AKANTU_NDEBUG)
+ int ret =
+#endif
+ MPI_Allreduce(MPI_IN_PLACE, values, nb_values, MPI_UNSIGNED,
+ synchronizer_operation_to_mpi_op[op],
+ communicator);
AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Allreduce.");
}
/* -------------------------------------------------------------------------- */
inline void StaticCommunicatorMPI::allReduce(Real * values, UInt nb_values,
const SynchronizerOperation & op) {
- int ret = MPI_Allreduce(MPI_IN_PLACE, values, nb_values, MPI_DOUBLE,
- synchronizer_operation_to_mpi_op[op],
- communicator);
+#if !defined(AKANTU_NDEBUG)
+ int ret =
+#endif
+ MPI_Allreduce(MPI_IN_PLACE, values, nb_values, MPI_DOUBLE,
+ synchronizer_operation_to_mpi_op[op],
+ communicator);
AKANTU_DEBUG_ASSERT(ret == MPI_SUCCESS, "Error in MPI_Allreduce.");
}
diff --git a/test/test_facet_extraction/test_facet_extraction_tetrahedron_4.cc b/test/test_facet_extraction/test_facet_extraction_tetrahedron_4.cc
index b59892712..5cf8a5de8 100644
--- a/test/test_facet_extraction/test_facet_extraction_tetrahedron_4.cc
+++ b/test/test_facet_extraction/test_facet_extraction_tetrahedron_4.cc
@@ -1,73 +1,73 @@
/**
* @file test_facet_extraction_tetra1.cc
* @author Guillaume ANCIAUX <anciaux@lsmscluster1.epfl.ch>
* @date Thu Aug 19 13:05:27 2010
*
* @brief test of internal facet extraction
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "mesh.hh"
#include "mesh_io.hh"
#include "mesh_io_msh.hh"
#include "mesh_utils.hh"
#include "solid_mechanics_model.hh"
#include "material.hh"
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_USE_IOHELPER
# include "io_helper.h"
#endif //AKANTU_USE_IOHELPER
using namespace akantu;
int main(int argc, char *argv[])
{
int dim = 3;
Mesh mesh(dim);
MeshIOMSH mesh_io;
mesh_io.read("cube.msh", mesh);
MeshUtils::buildFacets(mesh,1,1);
- unsigned int nb_nodes = mesh.getNbNodes();
#ifdef AKANTU_USE_IOHELPER
+ unsigned int nb_nodes = mesh.getNbNodes();
DumperParaview dumper;
dumper.SetMode(TEXT);
dumper.SetPoints(mesh.getNodes().values, dim, nb_nodes, "test-facet-extraction");
dumper.SetConnectivity((int*)mesh.getConnectivity(_tetrahedron_4).values,
TETRA1, mesh.getNbElement(_tetrahedron_4), C_MODE);
dumper.SetPrefix("paraview/");
dumper.Init();
dumper.Dump();
DumperParaview dumper_facet;
dumper_facet.SetMode(TEXT);
dumper_facet.SetPoints(mesh.getNodes().values, dim, nb_nodes, "test-facet-extraction_boundary");
dumper_facet.SetConnectivity((int*)mesh.getConnectivity(_triangle_3).values,
TRIANGLE1, mesh.getNbElement(_triangle_3), C_MODE);
dumper_facet.SetPrefix("paraview/");
dumper_facet.Init();
dumper_facet.Dump();
dumper_facet.SetPoints(mesh.getNodes().values, dim, nb_nodes, "test-facet-extraction_internal");
dumper_facet.SetConnectivity((int*)mesh.getInternalFacetsMesh().getConnectivity(_triangle_3).values,
TRIANGLE1, mesh.getInternalFacetsMesh().getNbElement(_triangle_3), C_MODE);
dumper_facet.Init();
dumper_facet.Dump();
#endif //AKANTU_USE_IOHELPER
return EXIT_SUCCESS;
}
diff --git a/test/test_facet_extraction/test_facet_extraction_triangle_3.cc b/test/test_facet_extraction/test_facet_extraction_triangle_3.cc
index 865a5ca24..e031d5782 100644
--- a/test/test_facet_extraction/test_facet_extraction_triangle_3.cc
+++ b/test/test_facet_extraction/test_facet_extraction_triangle_3.cc
@@ -1,72 +1,72 @@
/**
* @file test_facet_extraction.cc
* @author Guillaume ANCIAUX <anciaux@lsmscluster1.epfl.ch>
* @date Thu Aug 19 13:05:27 2010
*
* @brief test of internal facet extraction
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "mesh.hh"
#include "mesh_io.hh"
#include "mesh_io_msh.hh"
#include "mesh_utils.hh"
#include "solid_mechanics_model.hh"
#include "material.hh"
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_USE_IOHELPER
# include "io_helper.h"
#endif //AKANTU_USE_IOHELPER
using namespace akantu;
int main(int argc, char *argv[])
{
int dim = 2;
Mesh mesh(dim);
MeshIOMSH mesh_io;
mesh_io.read("square.msh", mesh);
MeshUtils::buildFacets(mesh,1,1);
- unsigned int nb_nodes = mesh.getNbNodes();
#ifdef AKANTU_USE_IOHELPER
+ unsigned int nb_nodes = mesh.getNbNodes();
DumperParaview dumper;
dumper.SetMode(TEXT);
dumper.SetPoints(mesh.getNodes().values, dim, nb_nodes, "test-facet-extraction");
dumper.SetConnectivity((int*)mesh.getConnectivity(_triangle_3).values,
TRIANGLE1, mesh.getNbElement(_triangle_3), C_MODE);
dumper.SetPrefix("paraview/");
dumper.Init();
dumper.Dump();
DumperParaview dumper_facet;
dumper_facet.SetMode(TEXT);
dumper_facet.SetPoints(mesh.getNodes().values, dim, nb_nodes, "test-facet-extraction_boundary");
dumper_facet.SetConnectivity((int*)mesh.getConnectivity(_segment_2).values,
LINE1, mesh.getNbElement(_segment_2), C_MODE);
dumper_facet.SetPrefix("paraview/");
dumper_facet.Init();
dumper_facet.Dump();
dumper_facet.SetPoints(mesh.getNodes().values, dim, nb_nodes, "test-facet-extraction_internal");
dumper_facet.SetConnectivity((int*)mesh.getInternalFacetsMesh().getConnectivity(_segment_2).values,
LINE1, mesh.getInternalFacetsMesh().getNbElement(_segment_2), C_MODE);
dumper_facet.Init();
dumper_facet.Dump();
#endif //AKANTU_USE_IOHELPER
return EXIT_SUCCESS;
}
diff --git a/test/test_mesh_io/test_mesh_io_msh.cc b/test/test_mesh_io/test_mesh_io_msh.cc
index 70d1f83bf..57886f9ad 100644
--- a/test/test_mesh_io/test_mesh_io_msh.cc
+++ b/test/test_mesh_io/test_mesh_io_msh.cc
@@ -1,34 +1,34 @@
/**
- * @file mesh_io_msh.cc
+ * @file test_mesh_io_msh.cc
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Wed Jul 14 14:27:11 2010
*
* @brief unit test for the MeshIOMSH class
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -------------------------------------------------------------------------- */
#include <cstdlib>
/* -------------------------------------------------------------------------- */
#include "mesh.hh"
#include "mesh_io.hh"
#include "mesh_io_msh.hh"
/* -------------------------------------------------------------------------- */
int main(int argc, char *argv[]) {
akantu::MeshIOMSH mesh_io;
akantu::Mesh mesh(3);
mesh_io.read("./cube.msh", mesh);
std::cout << mesh << std::endl;
return EXIT_SUCCESS;
}
diff --git a/test/test_mesh_partitionate/test_mesh_partitionate_scotch.cc b/test/test_mesh_partitionate/test_mesh_partitionate_scotch.cc
index 5b374345b..930819dd3 100644
--- a/test/test_mesh_partitionate/test_mesh_partitionate_scotch.cc
+++ b/test/test_mesh_partitionate/test_mesh_partitionate_scotch.cc
@@ -1,79 +1,79 @@
/**
- * @file test_mesh_partition_scotch.cc
+ * @file test_mesh_partitionate_scotch.cc
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Thu Aug 19 13:05:27 2010
*
* @brief test of internal facet extraction
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "mesh.hh"
#include "fem.hh"
#include "mesh_io_msh.hh"
#include "mesh_partition_scotch.hh"
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_USE_IOHELPER
# include "io_helper.h"
#endif //AKANTU_USE_IOHELPER
/* -------------------------------------------------------------------------- */
/* Main */
/* -------------------------------------------------------------------------- */
int main(int argc, char *argv[])
{
akantu::initialize(&argc, &argv);
akantu::debug::setDebugLevel(akantu::dblDump);
int dim = 2;
#ifdef AKANTU_USE_IOHELPER
akantu::ElementType type = akantu::_triangle_6;
#endif //AKANTU_USE_IOHELPER
akantu::Mesh mesh(dim);
akantu::MeshIOMSH mesh_io;
mesh_io.read("triangle.msh", mesh);
akantu::MeshPartition * partition = new akantu::MeshPartitionScotch(mesh, dim);
partition->partitionate(8);
#ifdef AKANTU_USE_IOHELPER
unsigned int nb_nodes = mesh.getNbNodes();
unsigned int nb_element = mesh.getNbElement(type);
DumperParaview dumper;
dumper.SetMode(TEXT);
dumper.SetPoints(mesh.getNodes().values, dim, nb_nodes, "test-scotch-partition");
dumper.SetConnectivity((int*) mesh.getConnectivity(type).values,
TRIANGLE2, nb_element, C_MODE);
akantu::UInt nb_quadrature_points = akantu::FEM::getNbQuadraturePoints(type);
double * part = new double[nb_element*nb_quadrature_points];
akantu::UInt * part_val = partition->getPartition(type).values;
for (unsigned int i = 0; i < nb_element; ++i)
for (unsigned int q = 0; q < nb_quadrature_points; ++q)
part[i*nb_quadrature_points + q] = part_val[i];
dumper.AddElemDataField(part, 1, "partitions");
dumper.SetPrefix("paraview/");
dumper.Init();
dumper.Dump();
delete [] part;
#endif //AKANTU_USE_IOHELPER
delete partition;
akantu::finalize();
return EXIT_SUCCESS;
}
diff --git a/test/test_static_memory/test_static_memory.cc b/test/test_static_memory/test_static_memory.cc
index 83b683d53..12ed97df7 100644
--- a/test/test_static_memory/test_static_memory.cc
+++ b/test/test_static_memory/test_static_memory.cc
@@ -1,32 +1,32 @@
/**
- * @file test/static_memory.cc
+ * @file test_static_memory.cc
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Fri Jun 11 11:55:54 2010
*
* @brief unit test for the StaticMemory class
*
*/
/* -------------------------------------------------------------------------- */
#include <iostream>
/* -------------------------------------------------------------------------- */
#include "aka_static_memory.hh"
#include "aka_vector.hh"
/* -------------------------------------------------------------------------- */
int main(int argc, char *argv[]) {
akantu::StaticMemory * st_mem = akantu::StaticMemory::getStaticMemory();
akantu::Vector<int> & test_int = st_mem->smalloc<int>(0, "test_int", 1000, 3);
test_int.resize(1050);
test_int.resize(2000);
std::cout << *st_mem << std::endl;
st_mem->sfree(0, "test_int");
exit(EXIT_SUCCESS);
}
diff --git a/test/test_vector/test_vector.cc b/test/test_vector/test_vector.cc
index d464d473e..a1b81e1fb 100644
--- a/test/test_vector/test_vector.cc
+++ b/test/test_vector/test_vector.cc
@@ -1,55 +1,55 @@
/**
- * @file vector.cc
+ * @file test_vector.cc
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Tue Jun 29 17:32:23 2010
*
* @brief test of the vector class
*
* @section LICENSE
*
* \<insert license here\>
*
*/
/* -------------------------------------------------------------------------- */
#include <cstdlib>
/* -------------------------------------------------------------------------- */
#include "aka_vector.hh"
/* -------------------------------------------------------------------------- */
int main(int argc, char *argv[]) {
int def_value[3];
def_value[0] = 10;
def_value[1] = 20;
def_value[2] = 30;
akantu::Vector<int> int_vect(10, 3, def_value, "test1");
for (unsigned int i = 5; i < int_vect.getSize(); ++i) {
for (unsigned int j = 0; j < int_vect.getNbComponent(); ++j) {
int_vect.values[i*int_vect.getNbComponent() + j] = def_value[j]*10;
}
}
std::cerr << int_vect;
int new_elem[3];
new_elem[0] = 1;
new_elem[1] = 2;
new_elem[2] = 3;
int_vect.push_back(new_elem);
int_vect.push_back(200);
int_vect.erase(0);
std::cerr << int_vect;
akantu::Vector<int> int_vect0(0, 3, def_value, "test2");
std::cerr << int_vect0;
int_vect0.push_back(new_elem);
std::cerr << int_vect0;
return EXIT_SUCCESS;
}