diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh
index e808fb0ab..653e2f1bb 100644
--- a/src/common/aka_common.hh
+++ b/src/common/aka_common.hh
@@ -1,412 +1,412 @@
/**
* @file aka_common.hh
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Fri Jun 11 09:48:06 2010
*
* @namespace akantu
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
* @section DESCRIPTION
*
* All common things to be included in the projects files
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_COMMON_HH__
#define __AKANTU_COMMON_HH__
/* -------------------------------------------------------------------------- */
#include <list>
#include <limits>
/* -------------------------------------------------------------------------- */
#define __BEGIN_AKANTU__ namespace akantu {
#define __END_AKANTU__ }
/* -------------------------------------------------------------------------- */
#include "aka_config.hh"
#include "aka_error.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
/* Common types */
/* -------------------------------------------------------------------------- */
typedef double Real;
typedef unsigned int UInt;
typedef unsigned long long UInt64;
typedef signed int Int;
typedef std::string ID;
static const Real UINT_INIT_VALUE = 0;
#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;
/* -------------------------------------------------------------------------- */
/* Mesh/FEM/Model types */
/* -------------------------------------------------------------------------- */
typedef UInt Surface;
typedef std::pair<Surface, Surface> SurfacePair;
typedef std::list< SurfacePair > SurfacePairList;
/* -------------------------------------------------------------------------- */
/// @boost sequence of element to loop on in global tasks
#define AKANTU_REGULAR_ELEMENT_TYPE \
(_not_defined) \
(_segment_2) \
(_segment_3) \
(_triangle_3) \
(_triangle_6) \
(_tetrahedron_4) \
(_tetrahedron_10) \
(_quadrangle_4) \
(_quadrangle_8) \
(_hexahedron_8) \
(_point) \
(_bernoulli_beam_2)
#define AKANTU_COHESIVE_ELEMENT_TYPE \
(_cohesive_2d_4) \
(_cohesive_2d_6)
#define AKANTU_ALL_ELEMENT_TYPE \
AKANTU_REGULAR_ELEMENT_TYPE AKANTU_COHESIVE_ELEMENT_TYPE
/// @enum ElementType type of element potentially contained in a Mesh
enum ElementType {
_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
_quadrangle_4, /// first order quadrangle
_quadrangle_8, /// second order quadrangle
_hexahedron_8, /// first order hexahedron
_point, /// point only for some algorithm to be generic like mesh partitioning
_bernoulli_beam_2, /// bernoulli beam 2D
_cohesive_2d_4, /// first order 2D cohesive
_cohesive_2d_6, /// second order 2D cohesive
_max_element_type
};
enum ElementKind {
_ek_not_defined,
_ek_regular,
_ek_cohesive
};
/// enum AnalysisMethod type of solving method used to solve the equation of motion
enum AnalysisMethod {
_static,
_implicit_dynamic,
_explicit_dynamic
};
/// myfunction(double * position, double * stress/force, double * normal, unsigned int material_id)
typedef void (*BoundaryFunction)(double *,double *, double*, unsigned int);
/// @enum BoundaryFunctionType type of function passed for boundary conditions
enum BoundaryFunctionType {
_bft_stress,
_bft_traction
};
/// @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,
_ct_rigid = 3
};
enum ContactSearchType {
_cst_not_defined = 0,
_cst_2d_expli = 1,
_cst_expli = 2
};
enum ContactNeighborStructureType {
_cnst_not_defined = 0,
_cnst_regular_grid = 1,
_cnst_2d_grid = 2
};
/* -------------------------------------------------------------------------- */
/* Friction */
/* -------------------------------------------------------------------------- */
typedef ID FrictionID;
/* -------------------------------------------------------------------------- */
/* Ghosts handling */
/* -------------------------------------------------------------------------- */
typedef ID SynchronizerID;
/// @enum CommunicatorType type of communication method to use
enum CommunicatorType {
_communicator_mpi,
_communicator_dummy
};
/// @enum SynchronizationTag type of synchronizations
enum SynchronizationTag {
/// SolidMechanicsModel tags
_gst_smm_mass, /// synchronization of the SolidMechanicsModel.mass
_gst_smm_for_strain, /// synchronization of the SolidMechanicsModel.current_position
_gst_smm_boundary, /// synchronization of the boundary, forces, velocities and displacement
_gst_smm_uv, /// synchronization of the nodal velocities and displacement
_gst_smm_res, /// synchronization of the nodal residual
_gst_smm_init_mat, /// synchronization of the data to initialize materials
/// HeatTransfer tags
_gst_htm_capacity, /// synchronization of the nodal heat capacity
_gst_htm_temperature, /// synchronization of the nodal temperature
_gst_htm_gradient_temperature, /// synchronization of the element gradient temperature
/// Material non local
_gst_mnl_damage, /// synchronization of data to average in non local material
/// 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_INCLUDE_INLINE_IMPL
/* -------------------------------------------------------------------------- */
#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_CONST(name, variable, type) \
inline const Vector<type> & \
get##name (const ::akantu::ElementType & el_type, \
const GhostType & ghost_type = _not_ghost) const { \
AKANTU_DEBUG_IN(); \
\
AKANTU_DEBUG_OUT(); \
return variable(el_type, ghost_type); \
}
#define AKANTU_GET_MACRO_BY_ELEMENT_TYPE(name, variable, type) \
inline Vector<type> & \
get##name (const ::akantu::ElementType & el_type, \
const GhostType & ghost_type = _not_ghost) { \
AKANTU_DEBUG_IN(); \
\
AKANTU_DEBUG_OUT(); \
return variable(el_type, ghost_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 _quadrangle_4 : stream << "quadrangle_4" ; break;
case _quadrangle_8 : stream << "quadrangle_8" ; break;
case _hexahedron_8 : stream << "hexahedron_8" ; break;
case _bernoulli_beam_2 : stream << "bernoulli_beam_2"; break;
case _cohesive_2d_4 : stream << "cohesive_2d_4" ; break;
case _cohesive_2d_6 : stream << "cohesive_2d_6" ; break;
case _not_defined : stream << "undefined" ; break;
case _max_element_type : stream << "ElementType(" << (int) type << ")"; break;
case _point : stream << "point"; break;
}
return stream;
}
/// standard output stream operator for GhostType
inline std::ostream & operator <<(std::ostream & stream, GhostType type)
{
switch(type)
{
case _not_ghost : stream << "not_ghost"; break;
case _ghost : stream << "ghost" ; break;
case _casper : stream << "Casper the friendly ghost"; 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);
+inline std::string to_lower(const std::string & str);
/* -------------------------------------------------------------------------- */
inline void trim(std::string & to_trim);
__END_AKANTU__
#include "aka_common_inline_impl.cc"
/* -------------------------------------------------------------------------- */
// BOOST PART: TOUCH ONLY IF YOU KNOW WHAT YOU ARE DOING
#include <boost/preprocessor.hpp>
#define AKANTU_EXCLUDE_ELEMENT_TYPE(type) \
AKANTU_DEBUG_ERROR("Type (" << type << ") not handled by this function")
#define AKANTU_BOOST_CASE_MACRO(r,macro,type) \
case type : { macro(type); break; }
#define AKANTU_BOOST_ELEMENT_SWITCH(macro1, list1, macro2, list2) \
do { \
switch(type) { \
BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_CASE_MACRO, macro1, list1) \
BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_CASE_MACRO, macro2, list2) \
case _max_element_type: { \
AKANTU_DEBUG_ERROR("Wrong type : " << type); \
break; \
} \
} \
} while(0)
#define AKANTU_BOOST_ALL_ELEMENT_SWITCH(macro1) \
do { \
switch(type) { \
BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_CASE_MACRO, macro1, AKANTU_ALL_ELEMENT_TYPE) \
case _max_element_type: { \
AKANTU_DEBUG_ERROR("Wrong type : " << type); \
break; \
} \
} \
} while(0)
#define AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(macro) \
AKANTU_BOOST_ELEMENT_SWITCH(macro, \
AKANTU_REGULAR_ELEMENT_TYPE, \
AKANTU_EXCLUDE_ELEMENT_TYPE, \
AKANTU_COHESIVE_ELEMENT_TYPE)
#define AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(macro) \
AKANTU_BOOST_ELEMENT_SWITCH(macro, \
AKANTU_COHESIVE_ELEMENT_TYPE, \
AKANTU_EXCLUDE_ELEMENT_TYPE, \
AKANTU_REGULAR_ELEMENT_TYPE)
#define AKANTU_BOOST_LIST_MACRO(r,macro,type) \
macro(type)
#define AKANTU_BOOST_ELEMENT_LIST(macro,list) \
BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_LIST_MACRO,macro,list)
#define AKANTU_BOOST_ALL_ELEMENT_LIST(macro) \
AKANTU_BOOST_ELEMENT_LIST(macro, AKANTU_ALL_ELEMENT_TYPE)
#define AKANTU_BOOST_REGULAR_ELEMENT_LIST(macro) \
AKANTU_BOOST_ELEMENT_LIST(macro, AKANTU_REGULAR_ELEMENT_TYPE)
#define AKANTU_BOOST_COHESIVE_ELEMENT_LIST(macro) \
AKANTU_BOOST_ELEMENT_LIST(macro, AKANTU_COHESIVE_ELEMENT_TYPE)
#endif /* __AKANTU_COMMON_HH__ */
diff --git a/src/common/aka_common_inline_impl.cc b/src/common/aka_common_inline_impl.cc
index c1432d8ba..6d660d413 100644
--- a/src/common/aka_common_inline_impl.cc
+++ b/src/common/aka_common_inline_impl.cc
@@ -1,54 +1,56 @@
/**
* @file aka_common_inline_impl.cc
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @date Fri Jun 11 09:48:06 2010
*
* @namespace akantu
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
* @section DESCRIPTION
*
* All common things to be included in the projects files
*
*/
#include <algorithm>
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
-inline void to_lower(std::string & str) {
+inline std::string to_lower(const std::string & str) {
+ std::string lstr;
std::transform(str.begin(),
str.end(),
- str.begin(),
+ lstr.begin(),
(int(*)(int))std::tolower);
+ return lstr;
}
/* -------------------------------------------------------------------------- */
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__
diff --git a/src/common/aka_error.cc b/src/common/aka_error.cc
index cc2187e90..fa7356a07 100644
--- a/src/common/aka_error.cc
+++ b/src/common/aka_error.cc
@@ -1,201 +1,201 @@
/**
* @file aka_error.cc
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Sun Sep 5 21:03:37 2010
*
* @brief
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_config.hh"
#include "aka_error.hh"
/* -------------------------------------------------------------------------- */
#include <iostream>
#include <csignal>
#include <execinfo.h>
#include <cxxabi.h>
#include <fstream>
#include <iomanip>
#include <cmath>
#include <sys/wait.h>
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
namespace debug {
/* ------------------------------------------------------------------------ */
void initSignalHandler() {
struct sigaction action;
action.sa_handler = &printBacktrace;
sigemptyset(&(action.sa_mask));
action.sa_flags = SA_RESETHAND;
sigaction(SIGSEGV, &action, NULL);
}
/* ------------------------------------------------------------------------ */
std::string demangle(const char* symbol) {
int status;
std::string result;
char * demangled_name;
if ((demangled_name = abi::__cxa_demangle(symbol, NULL, 0, &status)) != NULL) {
result = demangled_name;
free(demangled_name);
} else {
result = symbol;
}
return result;
//return symbol;
}
/* ------------------------------------------------------------------------ */
- void printBacktrace(int sig) {
+ void printBacktrace(__attribute__((unused)) int sig) {
AKANTU_DEBUG_INFO("Caught signal " << sig << "!");
// std::stringstream pidsstr;
// pidsstr << getpid();
// char name_buf[512];
// name_buf[readlink("/proc/self/exe", name_buf, 511)]=0;
// std::string execname(name_buf);
// std::cout << "stack trace for " << execname << " pid=" << pidsstr.str() << std::endl;
// std::string cmd;
// cmd = "CMDFILE=$(mktemp); echo 'bt' > ${CMDFILE}; gdb --batch " + execname + " " + pidsstr.str() + " < ${CMDFILE};";
// int retval __attribute__((unused)) = system(("bash -c '" + cmd + "'").c_str());
const size_t max_depth = 100;
size_t stack_depth;
void *stack_addrs[max_depth];
char **stack_strings;
size_t i;
stack_depth = backtrace(stack_addrs, max_depth);
stack_strings = backtrace_symbols(stack_addrs, stack_depth);
std::cerr << "BACKTRACE : " << stack_depth << " stack frames." <<std::endl;
size_t w = size_t(std::floor(log(double(stack_depth))/std::log(10.))+1);
/// -1 to remove the call to the printBacktrace function
for (i = 1; i < stack_depth; i++) {
std::cerr << " [" << std::setw(w) << i << "] ";
std::string bt_line(stack_strings[i]);
size_t first, second;
if((first = bt_line.find('(')) != std::string::npos && (second = bt_line.find('+')) != std::string::npos) {
std::cerr << bt_line.substr(0,first + 1) << demangle(bt_line.substr(first + 1, second - first - 1).c_str()) << bt_line.substr(second) << std::endl;
// char name_exe[512];
// name_exe[readlink("/proc/self/exe", name_exe, 511)]=0;
// std::stringstream syscom;
// syscom << "addr2line " << stack_addrs[i] << "-C -s -e " << name_exe;
// for (UInt i = 0; i < w + 4; ++i) std::cerr << " ";
// std::cout << "-> " << std::flush;
// int retval __attribute__((unused)) = system(syscom.str().c_str());
} else {
std::cerr << bt_line << std::endl;
}
}
free(stack_strings);
std::cerr << "END BACKTRACE" << std::endl;
}
/* ------------------------------------------------------------------------ */
/* ------------------------------------------------------------------------ */
Debugger::Debugger() {
cout = &std::cerr;
level = dblInfo;
parallel_context = "";
file_open = false;
}
/* ------------------------------------------------------------------------ */
Debugger::~Debugger() {
if(file_open) {
dynamic_cast<std::ofstream *>(cout)->close();
delete cout;
}
}
/* ------------------------------------------------------------------------ */
void Debugger::throwException(const std::string & info) throw(akantu::debug::Exception) {
AKANTU_DEBUG(akantu::dblWarning, "!!! " << info);
::akantu::debug::Exception ex(info, __FILE__, __LINE__ );
throw ex;
}
/* ------------------------------------------------------------------------ */
void Debugger::exit(int status) {
int * a = NULL;
*a = 1;
if (status != EXIT_SUCCESS)
akantu::debug::printBacktrace(15);
#ifdef AKANTU_USE_MPI
MPI_Abort(MPI_COMM_WORLD, MPI_ERR_UNKNOWN);
#endif
exit(status); // not called when compiled with MPI due to MPI_Abort, but
// MPI_Abort does not have the noreturn attribute
}
/* ------------------------------------------------------------------------ */
void Debugger::setDebugLevel(const DebugLevel & level) {
this->level = level;
}
/* ------------------------------------------------------------------------ */
const DebugLevel & Debugger::getDebugLevel() const {
return level;
}
/* ------------------------------------------------------------------------ */
void Debugger::setLogFile(const std::string & filename) {
if(file_open) {
dynamic_cast<std::ofstream *>(cout)->close();
delete cout;
}
std::ofstream * fileout = new std::ofstream(filename.c_str());
file_open = true;
cout = fileout;
}
std::ostream & Debugger::getOutputStream() {
return *cout;
}
/* ------------------------------------------------------------------------ */
void Debugger::setParallelContext(int rank, int size) {
std::stringstream sstr;
sstr << "[" << std::setfill(' ') << std::right << std::setw(3)
<< (rank + 1) << "/" << size << "] ";
parallel_context = sstr.str();
}
void setDebugLevel(const DebugLevel & level) {
debugger.setDebugLevel(level);
}
}
__END_AKANTU__
diff --git a/src/common/aka_error.hh b/src/common/aka_error.hh
index 80746571c..15ca59570 100644
--- a/src/common/aka_error.hh
+++ b/src/common/aka_error.hh
@@ -1,272 +1,280 @@
/**
* @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
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_ERROR_HH__
#define __AKANTU_ERROR_HH__
/* -------------------------------------------------------------------------- */
#include <ostream>
#include <sstream>
#include <sys/time.h>
#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,
dblDebug = 42,
dbl100 = 100,
dblDump = 100,
dblTest = 1337
};
/* -------------------------------------------------------------------------- */
namespace debug {
void setDebugLevel(const DebugLevel & level);
void initSignalHandler();
std::string demangle(const char * symbol);
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;
};
+ /// standard output stream operator
+ inline std::ostream & operator <<(std::ostream & stream, const Exception & _this)
+ {
+ stream << _this.what();
+ return stream;
+ }
+
+
/* -------------------------------------------------------------------------- */
#define AKANTU_LOCATION "(" <<__FILE__ << ":" << __func__ << "():" << __LINE__ << ")"
class Debugger {
public:
Debugger();
virtual ~Debugger();
void exit(int status) __attribute__ ((noreturn));
- void throwException(const std::string & info) throw(akantu::debug::Exception);
+ void throwException(const std::string & info) throw(akantu::debug::Exception) __attribute__ ((noreturn));
inline void printMessage(const std::string & prefix,
const DebugLevel & level,
const std::string & info) const {
if(this->level >= level) {
struct timeval time;
gettimeofday(&time, NULL);
double timestamp = time.tv_sec*1e6 + time.tv_usec; /*in us*/
*(cout) << parallel_context
<< "{" << (unsigned int)timestamp << "} "
<< prefix << info
<< std::endl;
}
}
public:
void setParallelContext(int rank, int size);
void setDebugLevel(const DebugLevel & level);
const DebugLevel & getDebugLevel() const;
void setLogFile(const std::string & filename);
std::ostream & getOutputStream();
inline bool testLevel(const DebugLevel & level) const {
return (this->level >= (level));
}
private:
std::string parallel_context;
std::ostream * cout;
bool file_open;
DebugLevel level;
};
extern Debugger debugger;
}
/* -------------------------------------------------------------------------- */
#define AKANTU_STRINGSTREAM_IN(_str, _sstr); \
do { \
std::stringstream _dbg_s_info; \
_dbg_s_info << _sstr; \
_str = _dbg_s_info.str(); \
} while(0)
/* -------------------------------------------------------------------------- */
#define AKANTU_EXCEPTION(info) \
do { \
std::string _dbg_str; \
AKANTU_STRINGSTREAM_IN(_dbg_str, info); \
::akantu::debug::debugger.throwException(_dbg_str); \
} while(0)
#define AKANTU_DEBUG_TO_IMPLEMENT() \
AKANTU_DEBUG_ERROR(__func__ << " : not implemented yet !"); \
exit(-1);
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_NDEBUG
#define AKANTU_DEBUG_TEST(level) (false)
#define AKANTU_DEBUG_LEVEL_IS_TEST() (false)
#define AKANTU_DEBUG(level,info)
#define AKANTU_DEBUG_(pref,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_(" ",level, info)
#define AKANTU_DEBUG_(pref, level, info) \
do { \
std::string _dbg_str; \
AKANTU_STRINGSTREAM_IN(_dbg_str, info << " " << AKANTU_LOCATION); \
::akantu::debug::debugger.printMessage(pref, level, _dbg_str); \
} while(0)
#define AKANTU_DEBUG_TEST(level) \
(::akantu::debug::debugger.testLevel(level))
#define AKANTU_DEBUG_LEVEL_IS_TEST() \
(::akantu::debug::debugger.testLevel(dblTest))
#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::debug::debugger.exit(EXIT_FAILURE); \
} \
} while(0)
#define AKANTU_DEBUG_ERROR(info) \
do { \
AKANTU_DEBUG_("!!! ", ::akantu::dblError, info); \
::akantu::debug::debugger.exit(EXIT_FAILURE); \
} while(0)
#endif // AKANTU_NDEBUG
/* -------------------------------------------------------------------------- */
__END_AKANTU__
#endif /* __AKANTU_ERROR_HH__ */
// LocalWords: acessory
diff --git a/src/common/aka_memory.hh b/src/common/aka_memory.hh
index 97b7e9aec..0d9b8a370 100644
--- a/src/common/aka_memory.hh
+++ b/src/common/aka_memory.hh
@@ -1,115 +1,118 @@
/**
* @file aka_memory.hh
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Tue Jun 15 09:30:23 2010
*
* @brief wrapper for the static_memory, all object which wants
* to access the ststic_memory as to inherit from the class memory
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MEMORY_HH__
#define __AKANTU_MEMORY_HH__
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "aka_static_memory.hh"
#include "aka_vector.hh"
/* -------------------------------------------------------------------------- */
#include <list>
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
class Memory {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
Memory(MemoryID memory_id = 0);
virtual ~Memory();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// malloc
template<class T>
inline Vector<T> & alloc(const ID & name,
UInt size,
UInt nb_component);
/// malloc
template<class T>
inline Vector<T> & alloc(const ID & name,
UInt size,
UInt nb_component,
const T & init_value);
/* ------------------------------------------------------------------------ */
/// free an array
inline void dealloc(const ID & name);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
AKANTU_GET_MACRO(MemoryID, memory_id, const MemoryID &);
- template<class T>
+ template<typename T>
inline Vector<T> & getVector(const ID & name);
+ template<typename T>
+ inline const Vector<T> & getVector(const ID & name) const;
+
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
/// the static memory instance
StaticMemory * static_memory;
protected:
/// the id registred in the static memory
MemoryID memory_id;
/// list of allocated vectors id
std::list<ID> handeld_vectors_id;
};
/* -------------------------------------------------------------------------- */
/* Inline functions */
/* -------------------------------------------------------------------------- */
#if defined (AKANTU_INCLUDE_INLINE_IMPL)
# include "aka_memory_inline_impl.cc"
#endif
__END_AKANTU__
#endif /* __AKANTU_MEMORY_HH__ */
diff --git a/src/common/aka_memory_inline_impl.cc b/src/common/aka_memory_inline_impl.cc
index d1600c577..afbac71a7 100644
--- a/src/common/aka_memory_inline_impl.cc
+++ b/src/common/aka_memory_inline_impl.cc
@@ -1,57 +1,62 @@
/**
* @file aka_memory_inline_impl.cc
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Thu Jul 15 00:22:48 2010
*
* @brief Implementation of the inline functions of the class Memory
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
template<class T> inline Vector<T> & Memory::alloc(const ID & name,
UInt size,
UInt nb_component) {
handeld_vectors_id.push_back(name);
return static_memory->smalloc<T>(memory_id, name,
size, nb_component);
}
/* -------------------------------------------------------------------------- */
template<class T> inline Vector<T> & Memory::alloc(const ID & name,
UInt size,
UInt nb_component,
const T & init_value) {
handeld_vectors_id.push_back(name);
return static_memory->smalloc<T>(memory_id, name,
size, nb_component, init_value);
}
/* -------------------------------------------------------------------------- */
inline void Memory::dealloc(const ID & name) {
AKANTU_DEBUG(dblAccessory, "Deleting the vector " << name);
static_memory->sfree(memory_id, name);
handeld_vectors_id.remove(name);
}
/* -------------------------------------------------------------------------- */
template<class T> inline Vector<T> & Memory::getVector(const ID & name) {
return static_cast< Vector<T> & >(const_cast<VectorBase &>(static_memory->getVector(memory_id, name)));
}
+
+/* -------------------------------------------------------------------------- */
+template<class T> inline const Vector<T> & Memory::getVector(const ID & name) const {
+ return static_cast< Vector<T> & >(const_cast<VectorBase &>(static_memory->getVector(memory_id, name)));
+}
diff --git a/src/fem/element_class_inline_impl.cc b/src/fem/element_class_inline_impl.cc
index 863172570..08d9a2861 100644
--- a/src/fem/element_class_inline_impl.cc
+++ b/src/fem/element_class_inline_impl.cc
@@ -1,409 +1,406 @@
/**
* @file element_class_inline_impl.cc
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @date Thu Jul 15 10:28:28 2010
*
* @brief Implementation of the inline functions of the class element_class
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
template<ElementType type> inline Real * ElementClass<type>::getQuadraturePoints() {
return quad;
}
/* -------------------------------------------------------------------------- */
template<ElementType type> inline UInt ElementClass<type>::getShapeSize() {
return nb_nodes_per_element;
}
/* -------------------------------------------------------------------------- */
template<ElementType type> inline UInt ElementClass<type>::getShapeDerivativesSize() {
return nb_nodes_per_element * spatial_dimension;
}
/* -------------------------------------------------------------------------- */
template<ElementType type>
void ElementClass<type>::preComputeStandards(const Real * coord,
const UInt dimension,
Real * shape,
Real * dshape,
Real * jacobians) {
// ask for computation of shapes
computeShapes(quad, nb_quadrature_points, shape);
// compute dnds
Real dnds[nb_nodes_per_element * spatial_dimension * nb_quadrature_points];
computeDNDS(quad, nb_quadrature_points, dnds);
// compute dxds
Real dxds[dimension * spatial_dimension * nb_quadrature_points];
computeDXDS(dnds, nb_quadrature_points, coord, dimension, dxds);
// jacobian
computeJacobian(dxds, nb_quadrature_points, dimension, jacobians);
// if dimension == spatial_dimension compute shape derivatives
if (dimension == spatial_dimension) {
computeShapeDerivatives(dxds, dnds, nb_quadrature_points, dimension, dshape);
}
}
/* -------------------------------------------------------------------------- */
template <ElementType type>
inline void ElementClass<type>::computeShapes(const Real * natural_coords,
const UInt nb_points,
Real * shapes) {
Real * cpoint = const_cast<Real *>(natural_coords);
for (UInt p = 0; p < nb_points; ++p) {
computeShapes(cpoint, shapes);
shapes += nb_nodes_per_element;
cpoint += spatial_dimension;
}
}
/* -------------------------------------------------------------------------- */
template <ElementType type>
inline void ElementClass<type>::computeShapes(const Real * natural_coords,
const UInt nb_points,
Real * shapes,
const Real * local_coord,
UInt id) {
Real * cpoint = const_cast<Real *>(natural_coords);
for (UInt p = 0; p < nb_points; ++p) {
computeShapes(cpoint, shapes, local_coord, id);
shapes += nb_nodes_per_element;
cpoint += spatial_dimension;
}
}
/* -------------------------------------------------------------------------- */
template<ElementType type>
inline void ElementClass<type>::computeDNDS(const Real * natural_coords,
const UInt nb_points,
Real * dnds) {
Real * cpoint = const_cast<Real *>(natural_coords);
Real * cdnds = dnds;
for (UInt p = 0; p < nb_points; ++p) {
computeDNDS(cpoint, cdnds);
cpoint += spatial_dimension;
cdnds += nb_nodes_per_element * spatial_dimension;
}
}
/* -------------------------------------------------------------------------- */
template<ElementType type>
inline void ElementClass<type>::computeDXDS(const Real * dnds,
const UInt nb_points,
const Real * node_coords,
const UInt dimension,
Real * dxds) {
Real * cdnds = const_cast<Real *>(dnds);
Real * cdxds = dxds;
for (UInt p = 0; p < nb_points; ++p) {
computeDXDS(cdnds, node_coords, dimension, cdxds);
cdnds += nb_nodes_per_element * spatial_dimension;
cdxds += spatial_dimension * dimension;
}
}
/* -------------------------------------------------------------------------- */
template <ElementType type>
inline void ElementClass<type>::computeDXDS(const Real * dnds,
const Real * node_coords,
const UInt dimension,
Real * dxds) {
/// @f$ J = dxds = dnds * x @f$
Math::matrix_matrix(spatial_dimension, dimension, nb_nodes_per_element,
dnds, node_coords, dxds);
}
/* -------------------------------------------------------------------------- */
template <ElementType type>
inline void ElementClass<type>::computeJacobian(const Real * dxds,
const UInt nb_points,
const UInt dimension,
Real * jac) {
Real * cdxds = const_cast<Real *>(dxds);
Real * cjac = jac;
for (UInt p = 0; p < nb_points; ++p) {
computeJacobian(cdxds, dimension, *cjac);
// AKANTU_DEBUG_ASSERT((cjac[0] > 0),
// "Negative jacobian computed, possible problem in the element node order.");
cdxds += spatial_dimension * dimension;
cjac++;
}
}
/* -------------------------------------------------------------------------- */
template <ElementType type>
inline void ElementClass<type>::computeShapeDerivatives(const Real * dxds,
const Real * dnds,
const UInt nb_points,
- const UInt dimension,
+ __attribute__ ((unused)) const UInt dimension,
Real * shape_deriv) {
AKANTU_DEBUG_ASSERT(dimension == spatial_dimension,"gradient in space "
<< dimension
<< " cannot be evaluated for element of dimension "
<< spatial_dimension);
Real * cdxds = const_cast<Real *>(dxds);
Real * cdnds = const_cast<Real *>(dnds);
for (UInt p = 0; p < nb_points; ++p) {
computeShapeDerivatives(cdxds, cdnds, shape_deriv);
cdnds += spatial_dimension * nb_nodes_per_element;
cdxds += spatial_dimension * spatial_dimension;
shape_deriv += nb_nodes_per_element * spatial_dimension;
}
}
/* -------------------------------------------------------------------------- */
template <ElementType type>
inline void ElementClass<type>::computeShapeDerivatives(const Real * natural_coords,
const UInt nb_points,
const UInt dimension,
Real * shape_deriv,
const Real * local_coord,
UInt id) {
AKANTU_DEBUG_ASSERT(dimension == spatial_dimension,"Gradient in space "
<< dimension
<< " cannot be evaluated for element of dimension "
<< spatial_dimension);
Real * cpoint = const_cast<Real *>(natural_coords);
for (UInt p = 0; p < nb_points; ++p) {
computeShapeDerivatives(cpoint, shape_deriv, local_coord, id);
shape_deriv += nb_nodes_per_element * dimension;
cpoint += dimension;
}
}
/* -------------------------------------------------------------------------- */
template <ElementType type>
inline void ElementClass<type>::computeShapeDerivatives(const Real * dxds,
const Real * dnds,
Real * shape_deriv) {
/// @f$ dxds = J^{-1} @f$
Real inv_dxds[spatial_dimension * spatial_dimension];
if (spatial_dimension == 1) inv_dxds[0] = 1./dxds[0];
if (spatial_dimension == 2) Math::inv2(dxds, inv_dxds);
if (spatial_dimension == 3) Math::inv3(dxds, inv_dxds);
Math::matrixt_matrixt(nb_nodes_per_element, spatial_dimension, spatial_dimension,
dnds, inv_dxds, shape_deriv);
}
/* -------------------------------------------------------------------------- */
template<ElementType type>
inline Real ElementClass<type>::getInradius(__attribute__ ((unused)) const Real * coord) {
AKANTU_DEBUG_ERROR("Function not implemented for type : " << type);
return 0;
}
/* -------------------------------------------------------------------------- */
template<ElementType type>
inline void ElementClass<type>::computeNormalsOnQuadPoint(const Real * coord,
const UInt dimension,
Real * normals) {
AKANTU_DEBUG_ASSERT((dimension - 1) == spatial_dimension,
"cannot extract a normal because of dimension mismatch "
<< dimension << " " << spatial_dimension);
Real * cpoint = const_cast<Real *>(quad);
Real * cnormals = normals;
Real dnds[spatial_dimension*nb_nodes_per_element];
Real dxds[spatial_dimension*dimension];
for (UInt p = 0; p < nb_quadrature_points; ++p) {
computeDNDS(cpoint,dnds);
computeDXDS(dnds,coord,dimension,dxds);
if (dimension == 2) {
Math::normal2(dxds,cnormals);
}
if (dimension == 3){
Math::normal3(dxds,dxds+dimension,cnormals);
}
cpoint += spatial_dimension;
cnormals += dimension;
}
}
/* -------------------------------------------------------------------------- */
template <ElementType type>
inline void ElementClass<type>::interpolateOnNaturalCoordinates(const Real * natural_coords,
const Real * nodal_values,
UInt dimension,
Real * interpolated){
Real shapes[nb_nodes_per_element];
computeShapes(natural_coords,shapes);
Math::matrix_matrix(1, dimension, nb_nodes_per_element,
shapes, nodal_values, interpolated);
}
/* -------------------------------------------------------------------------- */
template <ElementType type>
inline void ElementClass<type>::inverseMap(const types::RVector & real_coords,
const types::Matrix & node_coords,
UInt dimension,
types::RVector & natural_coords,
Real tolerance){
//matric copy of the real_coords
types::Matrix mreal_coords(real_coords.storage(),1,spatial_dimension);
//initial guess
types::Matrix natural_guess(1,dimension,0.);
// realspace coordinates provided by initial guess
types::Matrix physical_guess(1,dimension);
// objective function f = real_coords - physical_guess
types::Matrix f(1,dimension);
// dnds computed on the natural_guess
types::Matrix dnds(nb_nodes_per_element,spatial_dimension);
// dxds computed on the natural_guess
types::Matrix dxds(spatial_dimension,dimension);
// transposed dxds computed on the natural_guess
types::Matrix dxds_t(dimension,spatial_dimension);
// G = dxds * dxds_t
types::Matrix G(spatial_dimension,spatial_dimension);
// Ginv = G^{-1}
types::Matrix Ginv(spatial_dimension,spatial_dimension);
// J = Ginv * dxds
types::Matrix J(spatial_dimension,dimension);
// dxi = \xi_{k+1} - \xi in the iterative process
types::Matrix dxi(1,spatial_dimension);
/* --------------------------- */
// init before iteration loop
/* --------------------------- */
// do interpolation
interpolateOnNaturalCoordinates(natural_guess.storage(),node_coords.storage(),dimension,physical_guess.storage());
// compute initial objective function value f = real_coords - physical_guess
f = physical_guess;
f*= -1.;
f+= mreal_coords;
// compute initial error
Real inverse_map_error = f.norm();
/* --------------------------- */
// iteration loop
/* --------------------------- */
while(tolerance < inverse_map_error)
{
//compute dxds
computeDNDS(natural_guess.storage(), dnds.storage());
computeDXDS(dnds.storage(),node_coords.storage(),dimension,dxds.storage());
//compute G
dxds_t = dxds;
dxds_t.transpose();
G.mul<false,false>(dxds,dxds_t);
// inverse G
if (spatial_dimension == 1) Ginv[0] = 1./G[0];
else if (spatial_dimension == 2) Math::inv2(G.storage(),Ginv.storage());
else if (spatial_dimension == 3) Math::inv3(G.storage(),Ginv.storage());
//compute J
J.mul<false,false>(Ginv,dxds);
//compute increment
dxi.mul<false,false>(f,J);
//update our guess
natural_guess += dxi;
//interpolate
interpolateOnNaturalCoordinates(natural_guess.storage(),node_coords.storage(),dimension,physical_guess.storage());
// compute error
f = physical_guess;
f*= -1.;
f+= mreal_coords;
inverse_map_error = f.norm();
}
memcpy(natural_coords.storage(),natural_guess.storage(),sizeof(Real)*natural_coords.size());
}
/* -------------------------------------------------------------------------- */
template <ElementType type>
inline void ElementClass<type>::computeShapes(const Real * natural_coords,
Real * shapes) {
computeShapes(natural_coords, shapes, NULL,0);
}
/* -------------------------------------------------------------------------- */
template <ElementType type>
inline void ElementClass<type>::computeShapes(__attribute__ ((unused)) const Real * natural_coords,
__attribute__ ((unused)) Real * shapes,
__attribute__ ((unused)) const Real * local_coord,
__attribute__ ((unused)) UInt id) {
AKANTU_DEBUG_ERROR("Function not implemented for type : " << type);
}
/* -------------------------------------------------------------------------- */
template <ElementType type>
inline void ElementClass<type>::computeDNDS(__attribute__((unused)) const Real * natural_coords,
__attribute__((unused)) Real * dnds){
// computeDNDS(natural_coords, dnds);
AKANTU_DEBUG_TO_IMPLEMENT();
}
/* -------------------------------------------------------------------------- */
template <ElementType type>
inline void ElementClass<type>::computeShapeDerivatives(__attribute__ ((unused)) const Real * natural_coords,
__attribute__ ((unused)) Real * shape_deriv,
__attribute__ ((unused)) const Real * local_coord,
__attribute__ ((unused)) UInt id) {
AKANTU_DEBUG_ERROR("Function not implemented for type : " << type);
}
/* -------------------------------------------------------------------------- */
template <ElementType type>
inline void ElementClass<type>::computeJacobian(__attribute__ ((unused)) const Real * dxds,
__attribute__ ((unused)) const UInt dimension,
__attribute__ ((unused)) Real & jac) {
-//AKANTU_DEBUG_ERROR("Function not implemented for type : " << type);
+ //AKANTU_DEBUG_ERROR("Function not implemented for type : " << type);
}
/* -------------------------------------------------------------------------- */
template <ElementType type>
inline Real * ElementClass<type>::getGaussIntegrationWeights() {
return gauss_integration_weights;
}
/* -------------------------------------------------------------------------- */
template <ElementType type>
inline bool ElementClass<type>::contains(__attribute__ ((unused)) const types::RVector & natural_coords) {
AKANTU_DEBUG_ERROR("Function not implemented for type : " << type);
- return true;
+ return false;
}
/* -------------------------------------------------------------------------- */
-
-
-
#include "element_classes/element_class_segment_2_inline_impl.cc"
#include "element_classes/element_class_segment_3_inline_impl.cc"
#include "element_classes/element_class_triangle_3_inline_impl.cc"
#include "element_classes/element_class_triangle_6_inline_impl.cc"
#include "element_classes/element_class_tetrahedron_4_inline_impl.cc"
#include "element_classes/element_class_tetrahedron_10_inline_impl.cc"
#include "element_classes/element_class_quadrangle_4_inline_impl.cc"
#include "element_classes/element_class_quadrangle_8_inline_impl.cc"
#include "element_classes/element_class_hexahedron_8_inline_impl.cc"
#include "element_classes/element_class_bernoulli_beam_2_inline_impl.cc"
diff --git a/src/fem/fem_template_inline_impl.cc b/src/fem/fem_template_inline_impl.cc
index c4c94ec0d..a2dacd862 100644
--- a/src/fem/fem_template_inline_impl.cc
+++ b/src/fem/fem_template_inline_impl.cc
@@ -1,183 +1,166 @@
/**
* @file fem_template_inline_impl.hh
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @date Tue May 22 11:21:40 2012
*
* @brief
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
template <typename Integ, typename Shape>
inline void FEMTemplate<Integ,Shape>::inverseMap(const types::RVector & real_coords,
UInt element,
const ElementType & type,
types::RVector & natural_coords,
const GhostType & ghost_type) const{
AKANTU_DEBUG_IN();
- Mesh::type_iterator it = mesh->firstType(element_dimension, ghost_type);
- Mesh::type_iterator end = mesh->lastType(element_dimension, ghost_type);
- for(; it != end; ++it) {
- ElementType type = *it;
-
-
#define INVERSE_MAP(type) \
shape_functions.template inverseMap<type>(real_coords,element,natural_coords,ghost_type);
AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(INVERSE_MAP);
#undef INVERSE_MAP
-}
+
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <>
inline void FEMTemplate<IntegratorCohesive<IntegratorGauss>,ShapeCohesive<ShapeLagrange> >
-::inverseMap(const types::RVector & real_coords,
- UInt element,
- const ElementType & type,
- types::RVector & natural_coords,
- const GhostType & ghost_type) const{
-
+::inverseMap(__attribute__((unused)) const types::RVector & real_coords,
+ __attribute__((unused)) UInt element,
+ __attribute__((unused)) const ElementType & type,
+ __attribute__((unused)) types::RVector & natural_coords,
+ __attribute__((unused)) const GhostType & ghost_type) const{
AKANTU_DEBUG_TO_IMPLEMENT();
}
/* -------------------------------------------------------------------------- */
template <>
inline void FEMTemplate<IntegratorGauss,ShapeLinked >
-::inverseMap(const types::RVector & real_coords,
- UInt element,
- const ElementType & type,
- types::RVector & natural_coords,
- const GhostType & ghost_type) const{
+::inverseMap(__attribute__((unused)) const types::RVector & real_coords,
+ __attribute__((unused)) UInt element,
+ __attribute__((unused)) const ElementType & type,
+ __attribute__((unused)) types::RVector & natural_coords,
+ __attribute__((unused)) const GhostType & ghost_type) const{
AKANTU_DEBUG_TO_IMPLEMENT();
}
/* -------------------------------------------------------------------------- */
template <typename Integ, typename Shape>
inline bool FEMTemplate<Integ,Shape>::contains(const types::RVector & real_coords,
UInt element,
const ElementType & type,
const GhostType & ghost_type) const{
AKANTU_DEBUG_IN();
- Mesh::type_iterator it = mesh->firstType(element_dimension, ghost_type);
- Mesh::type_iterator end = mesh->lastType(element_dimension, ghost_type);
- for(; it != end; ++it) {
- ElementType type = *it;
-
+ bool contain = false;
#define CONTAINS(type) \
- return shape_functions.template contains<type>(real_coords,element,ghost_type);
+ contain = shape_functions.template contains<type>(real_coords,element,ghost_type);
- AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(CONTAINS);
+ AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(CONTAINS);
#undef CONTAINS
-}
+
AKANTU_DEBUG_OUT();
- return false;
+ return contain;
}
/* -------------------------------------------------------------------------- */
template <>
inline bool FEMTemplate<IntegratorCohesive<IntegratorGauss>,ShapeCohesive<ShapeLagrange> >
-::contains(const types::RVector & real_coords,
- UInt element,
- const ElementType & type,
- const GhostType & ghost_type) const{
+::contains(__attribute__((unused)) const types::RVector & real_coords,
+ __attribute__((unused)) UInt element,
+ __attribute__((unused)) const ElementType & type,
+ __attribute__((unused)) const GhostType & ghost_type) const{
AKANTU_DEBUG_TO_IMPLEMENT();
}
/* -------------------------------------------------------------------------- */
template <>
inline bool FEMTemplate<IntegratorGauss,ShapeLinked >
-::contains(const types::RVector & real_coords,
- UInt element,
- const ElementType & type,
- const GhostType & ghost_type) const{
+::contains(__attribute__((unused)) const types::RVector & real_coords,
+ __attribute__((unused)) UInt element,
+ __attribute__((unused)) const ElementType & type,
+ __attribute__((unused)) const GhostType & ghost_type) const{
AKANTU_DEBUG_TO_IMPLEMENT();
}
/* -------------------------------------------------------------------------- */
template <typename Integ, typename Shape>
inline void FEMTemplate<Integ,Shape>::computeShapes(const types::RVector & real_coords,
UInt element,
const ElementType & type,
types::RVector & shapes,
const GhostType & ghost_type) const{
AKANTU_DEBUG_IN();
- Mesh::type_iterator it = mesh->firstType(element_dimension, ghost_type);
- Mesh::type_iterator end = mesh->lastType(element_dimension, ghost_type);
- for(; it != end; ++it) {
- ElementType type = *it;
-
-
#define COMPUTE_SHAPES(type) \
shape_functions.template computeShapes<type>(real_coords,element,shapes,ghost_type);
AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(COMPUTE_SHAPES);
#undef COMPUTE_SHAPES
-}
+
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <>
inline void FEMTemplate<IntegratorCohesive<IntegratorGauss>,ShapeCohesive<ShapeLagrange> >
-::computeShapes(const types::RVector & real_coords,
- UInt element,
- const ElementType & type,
- types::RVector & shapes,
- const GhostType & ghost_type) const{
+::computeShapes(__attribute__((unused)) const types::RVector & real_coords,
+ __attribute__((unused)) UInt element,
+ __attribute__((unused)) const ElementType & type,
+ __attribute__((unused)) types::RVector & shapes,
+ __attribute__((unused)) const GhostType & ghost_type) const{
AKANTU_DEBUG_TO_IMPLEMENT();
}
/* -------------------------------------------------------------------------- */
template <>
inline void FEMTemplate<IntegratorGauss,ShapeLinked >
-::computeShapes(const types::RVector & real_coords,
- UInt element,
- const ElementType & type,
- types::RVector & shapes,
- const GhostType & ghost_type) const{
+::computeShapes(__attribute__((unused)) const types::RVector & real_coords,
+ __attribute__((unused)) UInt element,
+ __attribute__((unused)) const ElementType & type,
+ __attribute__((unused)) types::RVector & shapes,
+ __attribute__((unused)) const GhostType & ghost_type) const{
AKANTU_DEBUG_TO_IMPLEMENT();
}
diff --git a/src/fem/integrator_cohesive_inline_impl.cc b/src/fem/integrator_cohesive_inline_impl.cc
index a2c7d3721..c22364422 100644
--- a/src/fem/integrator_cohesive_inline_impl.cc
+++ b/src/fem/integrator_cohesive_inline_impl.cc
@@ -1,403 +1,403 @@
/**
* @file integrator_cohesive_inline_impl.cc
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Wed Feb 22 17:22:21 2012
*
* @brief
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
template <class Inte>
IntegratorCohesive<Inte>::IntegratorCohesive(const Mesh & mesh,
const ID & id,
const MemoryID & memory_id) :
Inte(mesh, id, memory_id) {
AKANTU_DEBUG_IN();
std::stringstream sstr;
sstr << id << "sub_integrator";
//sub_type_integrator = new Inte(mesh, sstr.str(), memory_id);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<class Inte>
template <ElementType type>
void IntegratorCohesive<Inte>::precomputeJacobiansOnQuadraturePoints(const GhostType & ghost_type) {
AKANTU_DEBUG_IN();
const ElementType sub_type = ElementType(CohesiveElementSubElementType<type>::value);
// sub_type_integrator->precomputeJacobiansOnQuadraturePoints<sub_type>(ghost_type);
/* ---------------------------------*/
UInt spatial_dimension = Inte::mesh->getSpatialDimension();
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
UInt nb_quadrature_points = CohesiveElement<type>::getNbQuadraturePoints();
Real * weights = CohesiveElement<type>::getGaussIntegrationWeights();
UInt * elem_val = Inte::mesh->getConnectivity(type,ghost_type).storage();;
UInt nb_element = Inte::mesh->getConnectivity(type,ghost_type).getSize();
Vector<Real> & jacobians_tmp = Inte::jacobians.alloc(nb_element*nb_quadrature_points,
1,
type,
ghost_type);
Real * jacobians_val = jacobians_tmp.storage();
Real local_coord[spatial_dimension * nb_nodes_per_element/2];
for (UInt elem = 0; elem < nb_element; ++elem) {
// extract the coordinates of the first line nodes first
UInt * connectivity = elem_val+elem*nb_nodes_per_element;
for (UInt n = 0; n < nb_nodes_per_element/2; ++n) {
memcpy(local_coord + n * spatial_dimension,
Inte::mesh->getNodes().storage() + connectivity[n] * spatial_dimension,
spatial_dimension * sizeof(Real));
}
static_cast<Inte*>(this)->computeJacobianOnQuadPointsByElement<sub_type>(spatial_dimension,
local_coord,
nb_nodes_per_element/2,
jacobians_val);
for (UInt q = 0; q < nb_quadrature_points; ++q) {
*jacobians_val++ *= weights[q];
}
// jacobians_val += nb_quadrature_points;
}
/* ---------------------------------*/
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<class Inte>
template <ElementType type>
inline void IntegratorCohesive<Inte>::integrateOnElement(const Vector<Real> & f,
Real * intf,
UInt nb_degree_of_freedom,
const UInt elem,
const GhostType & ghost_type) const {
AKANTU_DEBUG_IN();
const ElementType sub_type = ElementType(CohesiveElementSubElementType<type>::value);
static_cast<Inte*>(this)->integrateOnElement<sub_type>(f, intf, nb_degree_of_freedom, elem, ghost_type);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<class Inte>
template <ElementType type>
const Vector<Real> &
IntegratorCohesive<Inte>::getQuadraturePoints(const GhostType & ghost_type) const {
AKANTU_DEBUG_IN();
const ElementType sub_type = ElementType(CohesiveElementSubElementType<type>::value);
AKANTU_DEBUG_OUT();
return Inte::template getQuadraturePoints<sub_type>(ghost_type);
}
/* -------------------------------------------------------------------------- */
template<class Inte>
template <ElementType type>
void IntegratorCohesive<Inte>::computeQuadraturePoints(const GhostType & ghost_type) {
AKANTU_DEBUG_IN();
const ElementType sub_type = ElementType(CohesiveElementSubElementType<type>::value);
static_cast<Inte*>(this)->template computeQuadraturePoints<sub_type>(ghost_type);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<class Inte>
template <ElementType type>
inline void IntegratorCohesive<Inte>::
computeJacobianOnQuadPointsByElement(UInt spatial_dimension,
Real * node_coords,
UInt nb_nodes_per_element,
Real * jacobians) {
AKANTU_DEBUG_IN();
const ElementType sub_type = ElementType(CohesiveElementSubElementType<type>::value);
static_cast<Inte*>(this)->computeJacobianOnQuadPointsByElement<sub_type>(spatial_dimension,
node_coords,
nb_nodes_per_element,
jacobians);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<class Inte>
template <ElementType type>
void IntegratorCohesive<Inte>::integrate(const Vector<Real> & in_f,
Vector<Real> &intf,
UInt nb_degree_of_freedom,
const GhostType & ghost_type,
const Vector<UInt> * filter_elements) const {
AKANTU_DEBUG_IN();
// const ElementType sub_type = ElementType(CohesiveElementSubElementType<type>::value);
// Inte::template integrate<sub_type>(in_f, intf, nb_degree_of_freedom, ghost_type, filter_elements);
/* ------------------------------------------------------------------------ */
AKANTU_DEBUG_ASSERT(Inte::jacobians.exists(type, ghost_type),
"No jacobians for the type "
<< Inte::jacobians.printType(type, ghost_type));
UInt nb_element = Inte::mesh->getNbElement(type,ghost_type);
const Vector<Real> & jac_loc = Inte::jacobians(type, ghost_type);
UInt nb_quadrature_points = CohesiveElement<type>::getNbQuadraturePoints();
UInt * filter_elem_val = NULL;
if(filter_elements != NULL) {
nb_element = filter_elements->getSize();
filter_elem_val = filter_elements->values;
}
Real * in_f_val = in_f.storage();
Real * intf_val = intf.storage();
Real * jac_val = jac_loc.storage();
UInt offset_in_f = in_f.getNbComponent()*nb_quadrature_points;
UInt offset_intf = intf.getNbComponent();
Real * jac = jac_val;
for (UInt el = 0; el < nb_element; ++el) {
if(filter_elements != NULL) {
jac = jac_val + filter_elem_val[el] * nb_quadrature_points;
}
Inte::integrate(in_f_val, jac, intf_val, nb_degree_of_freedom, nb_quadrature_points);
in_f_val += offset_in_f;
intf_val += offset_intf;
if(filter_elements == NULL) {
jac += nb_quadrature_points;
}
}
/* ------------------------------------------------------------------------ */
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<class Inte>
template <ElementType type>
Real IntegratorCohesive<Inte>::integrate(const Vector<Real> & in_f,
const GhostType & ghost_type,
const Vector<UInt> * filter_elements) const {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(Inte::jacobians.exists(type, ghost_type),
"No jacobians for the type "
<< Inte::jacobians.printType(type, ghost_type));
UInt nb_element = Inte::mesh->getNbElement(type, ghost_type);
const Vector<Real> & jac_loc = Inte::jacobians(type, ghost_type);
UInt nb_quadrature_points = CohesiveElement<type>::getNbQuadraturePoints();
UInt * filter_elem_val = NULL;
if(filter_elements != NULL) {
nb_element = filter_elements->getSize();
filter_elem_val = filter_elements->values;
}
Real intf = 0.;
Real * in_f_val = in_f.storage();
Real * jac_val = jac_loc.storage();
UInt offset_in_f = in_f.getNbComponent() * nb_quadrature_points;
Real * jac = jac_val;
for (UInt el = 0; el < nb_element; ++el) {
if(filter_elements != NULL) {
jac = jac_val + filter_elem_val[el] * nb_quadrature_points;
}
Real el_intf = 0;
Inte::integrate(in_f_val, jac, &el_intf, 1, nb_quadrature_points);
intf += el_intf;
in_f_val += offset_in_f;
if(filter_elements == NULL) {
jac += nb_quadrature_points;
}
}
AKANTU_DEBUG_OUT();
return intf;
}
/* -------------------------------------------------------------------------- */
// template<class Inte>
// template <ElementType type>
// Real IntegratorCohesive<Inte>::integrate(const Vector<Real> & in_f,
// const GhostType & ghost_type,
// const Vector<UInt> * filter_elements) const {
// AKANTU_DEBUG_IN();
// const ElementType sub_type = ElementType(CohesiveElementSubElementType<type>::value);
// AKANTU_DEBUG_OUT();
// return Inte::template integrate<sub_type>(in_f, ghost_type, filter_elements);
// }
/* -------------------------------------------------------------------------- */
template<class Inte>
template <ElementType type>
void IntegratorCohesive<Inte>::integrateOnQuadraturePoints(const Vector<Real> & in_f,
Vector<Real> &intf,
UInt nb_degree_of_freedom,
const GhostType & ghost_type,
const Vector<UInt> * filter_elements) const {
AKANTU_DEBUG_IN();
const ElementType sub_type = ElementType(CohesiveElementSubElementType<type>::value);
Inte::template integrateOnQuadraturePoints<sub_type>(in_f,
intf,
nb_degree_of_freedom,
ghost_type,
filter_elements);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<class Inte>
template<ElementType type>
void IntegratorCohesive<Inte>::checkJacobians(const GhostType & ghost_type) const {
AKANTU_DEBUG_IN();
// const ElementType sub_type = ElementType(CohesiveElementSubElementType<type>::value);
// sub_type_integrator->checkJacobians<sub_type>(ghost_type);
UInt nb_quadrature_points = CohesiveElement<type>::getNbQuadraturePoints();
UInt nb_element;
nb_element = Inte::mesh->getConnectivity(type,ghost_type).getSize();
Real * jacobians_val = Inte::jacobians(type, ghost_type).storage();
for (UInt i = 0; i < nb_element*nb_quadrature_points; ++i,++jacobians_val){
AKANTU_DEBUG_ASSERT(*jacobians_val >0,
"Negative jacobian computed,"
<< " possible problem in the element node order");
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/* specialization */
/* -------------------------------------------------------------------------- */
template<>
template<>
inline void IntegratorCohesive<IntegratorGauss>
-::precomputeJacobiansOnQuadraturePoints<_not_defined>(const GhostType & ghost_type) {
+::precomputeJacobiansOnQuadraturePoints<_not_defined>(__attribute__((unused)) const GhostType & ghost_type) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/* -------------------------------------------------------------------------- */
template<>
template<>
inline void IntegratorCohesive<IntegratorGauss>
-::integrateOnElement<_not_defined>(const Vector<Real> & f,
- Real * intf,
- UInt nb_degree_of_freedom,
- const UInt elem,
- const GhostType & ghost_type) const {
+::integrateOnElement<_not_defined>(__attribute__((unused)) const Vector<Real> & f,
+ __attribute__((unused)) Real * intf,
+ __attribute__((unused)) UInt nb_degree_of_freedom,
+ __attribute__((unused)) const UInt elem,
+ __attribute__((unused)) const GhostType & ghost_type) const {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/* -------------------------------------------------------------------------- */
template<>
template<>
inline const Vector<Real> &
IntegratorCohesive<IntegratorGauss>
-::getQuadraturePoints<_not_defined>(const GhostType & ghost_type) const {
+::getQuadraturePoints<_not_defined>(__attribute__((unused)) const GhostType & ghost_type) const {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/* -------------------------------------------------------------------------- */
template<>
template<>
inline void IntegratorCohesive<IntegratorGauss>
-::computeQuadraturePoints<_not_defined>(const GhostType & ghost_type) {
+::computeQuadraturePoints<_not_defined>(__attribute__((unused)) const GhostType & ghost_type) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/* -------------------------------------------------------------------------- */
template<>
template<>
inline void IntegratorCohesive<IntegratorGauss>::
computeJacobianOnQuadPointsByElement<_not_defined>(__attribute__((unused)) UInt spatial_dimension,
__attribute__((unused)) Real * node_coords,
__attribute__((unused)) UInt nb_nodes_per_element,
__attribute__((unused)) Real * jacobians) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/* -------------------------------------------------------------------------- */
template<>
template<>
inline void IntegratorCohesive<IntegratorGauss>
-::integrate<_not_defined>(const Vector<Real> & in_f,
- Vector<Real> &intf,
- UInt nb_degree_of_freedom,
- const GhostType & ghost_type,
- const Vector<UInt> * filter_elements) const {
+::integrate<_not_defined>(__attribute__((unused)) const Vector<Real> & in_f,
+ __attribute__((unused)) Vector<Real> &intf,
+ __attribute__((unused)) UInt nb_degree_of_freedom,
+ __attribute__((unused)) const GhostType & ghost_type,
+ __attribute__((unused)) const Vector<UInt> * filter_elements) const {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/* -------------------------------------------------------------------------- */
template<>
template<>
inline Real IntegratorCohesive<IntegratorGauss>
-::integrate<_not_defined>(const Vector<Real> & in_f,
- const GhostType & ghost_type,
- const Vector<UInt> * filter_elements) const {
+::integrate<_not_defined>(__attribute__((unused)) const Vector<Real> & in_f,
+ __attribute__((unused)) const GhostType & ghost_type,
+ __attribute__((unused)) const Vector<UInt> * filter_elements) const {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/* -------------------------------------------------------------------------- */
template<>
template<>
inline void IntegratorCohesive<IntegratorGauss>
-::integrateOnQuadraturePoints<_not_defined>(const Vector<Real> & in_f,
- Vector<Real> &intf,
- UInt nb_degree_of_freedom,
- const GhostType & ghost_type,
- const Vector<UInt> * filter_elements) const {
+::integrateOnQuadraturePoints<_not_defined>(__attribute__((unused)) const Vector<Real> & in_f,
+ __attribute__((unused)) Vector<Real> &intf,
+ __attribute__((unused)) UInt nb_degree_of_freedom,
+ __attribute__((unused)) const GhostType & ghost_type,
+ __attribute__((unused)) const Vector<UInt> * filter_elements) const {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/* -------------------------------------------------------------------------- */
template<>
template<>
inline void IntegratorCohesive<IntegratorGauss>
-::checkJacobians<_not_defined>(const GhostType & ghost_type) const {
+::checkJacobians<_not_defined>(__attribute__((unused)) const GhostType & ghost_type) const {
AKANTU_DEBUG_TO_IMPLEMENT();
}
diff --git a/src/model/model.cc b/src/model/model.cc
index cedd43519..03edf5f34 100644
--- a/src/model/model.cc
+++ b/src/model/model.cc
@@ -1,144 +1,148 @@
/**
* @file model.cc
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Fri Sep 16 14:46:01 2011
*
* @brief
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "model.hh"
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
Model::Model(const ID & id,
const MemoryID & memory_id) :
Memory(memory_id), id(id),
synch_registry(NULL), is_pbc_slave_node(0, 1, "is_pbc_slave_node") {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void Model::createSynchronizerRegistry(DataAccessor * data_accessor){
synch_registry = new SynchronizerRegistry(*data_accessor);
}
/* -------------------------------------------------------------------------- */
void Model::setPBC(UInt x, UInt y, UInt z){
Mesh & mesh = getFEM().getMesh();
mesh.computeBoundingBox();
if (x) MeshUtils::computePBCMap(mesh, 0, pbc_pair);
if (y) MeshUtils::computePBCMap(mesh, 1, pbc_pair);
if (z) MeshUtils::computePBCMap(mesh, 2, pbc_pair);
}
/* -------------------------------------------------------------------------- */
void Model::setPBC(SurfacePairList & surface_pairs,
ElementType surface_e_type){
Mesh & mesh = getFEM().getMesh();
SurfacePairList::iterator s_it;
for(s_it = surface_pairs.begin(); s_it != surface_pairs.end(); ++s_it) {
MeshUtils::computePBCMap(mesh, *s_it, surface_e_type, pbc_pair);
}
}
/* -------------------------------------------------------------------------- */
void Model::initPBC() {
Mesh & mesh = getFEM().getMesh();
std::map<UInt,UInt>::iterator it = pbc_pair.begin();
std::map<UInt,UInt>::iterator end = pbc_pair.end();
is_pbc_slave_node.resize(mesh.getNbNodes());
-
+#ifndef AKANTU_NDEBUG
Real * coords = mesh.getNodes().values;
UInt dim = mesh.getSpatialDimension();
+#endif
while(it != end){
UInt i1 = (*it).first;
+
+#ifndef AKANTU_NDEBUG
UInt i2 = (*it).second;
+#endif
is_pbc_slave_node(i1) = true;
AKANTU_DEBUG_INFO("pairing " << i1 << " ("
<< coords[dim*i1] << "," << coords[dim*i1+1] << ","
<< coords[dim*i1+2]
<< ") with "
<< i2 << " ("
<< coords[dim*i2] << "," << coords[dim*i2+1] << ","
<< coords[dim*i2+2]
<< ")");
++it;
}
}
/* -------------------------------------------------------------------------- */
Synchronizer & Model::createParallelSynch(MeshPartition * partition,
__attribute__((unused)) DataAccessor * data_accessor){
AKANTU_DEBUG_IN();
/* ------------------------------------------------------------------------ */
/* Parallel initialization */
/* ------------------------------------------------------------------------ */
StaticCommunicator * comm =
StaticCommunicator::getStaticCommunicator();
Int prank = comm->whoAmI();
DistributedSynchronizer * synch = NULL;
if(prank == 0)
synch =
DistributedSynchronizer::createDistributedSynchronizerMesh(getFEM().getMesh(),
partition);
else
synch =
DistributedSynchronizer::createDistributedSynchronizerMesh(getFEM().getMesh(),
NULL);
AKANTU_DEBUG_OUT();
return *synch;
}
/* -------------------------------------------------------------------------- */
Model::~Model() {
AKANTU_DEBUG_IN();
FEMMap::iterator it;
for (it = fems.begin(); it != fems.end(); ++it) {
if(it->second) delete it->second;
}
for (it = fems_boundary.begin(); it != fems_boundary.end(); ++it) {
if(it->second) delete it->second;
}
delete synch_registry;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
__END_AKANTU__
diff --git a/src/model/solid_mechanics/material.cc b/src/model/solid_mechanics/material.cc
index 413282243..094affb0d 100644
--- a/src/model/solid_mechanics/material.cc
+++ b/src/model/solid_mechanics/material.cc
@@ -1,536 +1,570 @@
/**
* @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
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "material.hh"
#include "solid_mechanics_model.hh"
#include "sparse_matrix.hh"
#include "dof_synchronizer.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
Material::Material(SolidMechanicsModel & model, const ID & id) :
Memory(model.getMemoryID()),
id(id),
name(""),
model(&model),
stress("stress", id),
strain("strain", id),
element_filter("element_filter", id),
// potential_energy_vector(false),
potential_energy("potential_energy", id),
is_non_local(false), is_init(false) {
AKANTU_DEBUG_IN();
rho = 0;
AKANTU_DEBUG_ASSERT(this->model,"model has wrong type: cannot proceed");
spatial_dimension = this->model->getSpatialDimension();
/// allocate strain stress for local elements
initInternalVector(strain, spatial_dimension * spatial_dimension);
initInternalVector(stress, spatial_dimension * spatial_dimension);
/// for each connectivity types allocate the element filer array of the material
initInternalVector(element_filter, 1);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
Material::~Material() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
bool Material::setParam(const std::string & key, const std::string & value,
__attribute__ ((unused)) const ID & id) {
std::stringstream sstr(value);
if(key == "name") name = std::string(value);
else if(key == "rho") { sstr >> rho; }
else return false;
return true;
}
/* -------------------------------------------------------------------------- */
void Material::initMaterial() {
AKANTU_DEBUG_IN();
resizeInternalVector(stress);
resizeInternalVector(strain);
is_init = true;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<typename T>
void Material::initInternalVector(ByElementTypeVector<T> & vect,
UInt nb_component,
ElementKind element_kind) const {
AKANTU_DEBUG_IN();
model->getFEM().getMesh().initByElementTypeVector(vect,
nb_component,
spatial_dimension,
false,
element_kind);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<typename T>
void Material::resizeInternalVector(ByElementTypeVector<T> & by_el_type_vect,
ElementKind element_kind) const {
AKANTU_DEBUG_IN();
FEM * fem = & model->getFEM();
if (element_kind == _ek_cohesive)
fem = & model->getFEM("CohesiveFEM");
const Mesh & mesh = fem->getMesh();
for(UInt g = _not_ghost; g <= _ghost; ++g) {
GhostType gt = (GhostType) g;
Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, element_kind);
Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, element_kind);
for(; it != end; ++it) {
const Vector<UInt> & elem_filter = element_filter(*it, gt);
UInt nb_element = elem_filter.getSize();
UInt nb_quadrature_points = fem->getNbQuadraturePoints(*it, gt);
UInt new_size = nb_element * nb_quadrature_points;
Vector<T> & vect = by_el_type_vect(*it, gt);
UInt size = vect.getSize();
UInt nb_component = vect.getNbComponent();
vect.resize(new_size);
memset(vect.storage() + size * nb_component, 0, (new_size - size) * nb_component * sizeof(T));
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Compute the residual by assembling @f$\int_{e} \sigma_e \frac{\partial
* \varphi}{\partial X} dX @f$
*
* @param[in] displacements nodes displacements
* @param[in] ghost_type compute the residual for _ghost or _not_ghost element
*/
void Material::updateResidual(Vector<Real> & displacement, GhostType ghost_type) {
AKANTU_DEBUG_IN();
computeStress(displacement, ghost_type);
assembleResidual(ghost_type);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void Material::assembleResidual(GhostType ghost_type) {
AKANTU_DEBUG_IN();
UInt spatial_dimension = model->getSpatialDimension();
Vector<Real> & residual = const_cast<Vector<Real> &>(model->getResidual());
Mesh & mesh = model->getFEM().getMesh();
Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type);
Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type);
for(; it != last_type; ++it) {
const Vector<Real> & shapes_derivatives = model->getFEM().getShapesDerivatives(*it, ghost_type);
Vector<UInt> & elem_filter = element_filter(*it, ghost_type);
UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent();
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it);
UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(*it, ghost_type);
UInt nb_element = elem_filter.getSize();
/// 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.storage();
Real * shapesd_val;
UInt * elem_filter_val = elem_filter.storage();
Vector<Real> * shapesd_filtered =
new Vector<Real>(nb_element*nb_quadrature_points, size_of_shapes_derivatives, "filtered shapesd");
Real * shapesd_filtered_val = shapesd_filtered->values;
for (UInt el = 0; el < nb_element; ++el) {
shapesd_val = shapesd + elem_filter_val[el] * size_of_shapes_derivatives * nb_quadrature_points;
memcpy(shapesd_filtered_val, shapesd_val,
size_of_shapes_derivatives * nb_quadrature_points * sizeof(Real));
shapesd_filtered_val += size_of_shapes_derivatives * nb_quadrature_points;
}
Vector<Real> & stress_vect = stress(*it, ghost_type);
// Vector<Real>::iterator<types::Matrix> sigma = stress_vect.begin(spatial_dimension, spatial_dimension);
// Vector<Real>::iterator<types::Matrix> sigma_end = stress_vect.end(spatial_dimension, spatial_dimension);
// Vector<Real>::iterator<types::Matrix> nabla_B = shapesd_filtered->begin(nb_nodes_per_element, spatial_dimension);
// Vector<Real>::iterator<types::Matrix> sigma_dphi_dx_it = sigma_dphi_dx->begin(nb_nodes_per_element, spatial_dimension);
// for (; sigma != sigma_end; ++sigma, ++nabla_B, ++sigma_dphi_dx_it) {
// sigma_dphi_dx_it->mul<true,false>(*nabla_B, *sigma);
// }
Math::matrix_matrixt(nb_nodes_per_element, spatial_dimension, spatial_dimension,
*shapesd_filtered,
stress_vect,
*sigma_dphi_dx);
delete shapesd_filtered;
/**
* 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>(nb_element, 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,
model->getDOFSynchronizer().getLocalDOFEquationNumbers(),
residual.getNbComponent(),
*it, ghost_type, &elem_filter, -1);
delete int_sigma_dphi_dx;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Compute the stress from the strain
*
* @param[in] current_position nodes postition + displacements
* @param[in] ghost_type compute the residual for _ghost or _not_ghost element
*/
void Material::computeStress(Vector<Real> & displacement, GhostType ghost_type) {
AKANTU_DEBUG_IN();
resizeInternalVector(stress);
resizeInternalVector(strain);
UInt spatial_dimension = model->getSpatialDimension();
Mesh::type_iterator it = model->getFEM().getMesh().firstType(spatial_dimension, ghost_type);
Mesh::type_iterator last_type = model->getFEM().getMesh().lastType(spatial_dimension, ghost_type);
for(; it != last_type; ++it) {
Vector<UInt> & elem_filter = element_filter(*it, ghost_type);
Vector<Real> & strain_vect = strain(*it, ghost_type);
UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(*it, ghost_type);
UInt nb_element = elem_filter.getSize();
strain_vect.resize(nb_quadrature_points * nb_element);
/// compute @f$\nabla u@f$
model->getFEM().gradientOnQuadraturePoints(displacement, strain_vect,
spatial_dimension,
*it, ghost_type, &elem_filter);
/// compute @f$\mathbf{\sigma}_q@f$ from @f$\nabla u@f$
computeStress(*it, ghost_type);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Compute the stiffness matrix by assembling @f$\int_{\omega} B^t \times D
* \times B d\omega @f$
*
* @param[in] current_position nodes postition + displacements
* @param[in] ghost_type compute the residual for _ghost or _not_ghost element
*/
void Material::assembleStiffnessMatrix(Vector<Real> & current_position, GhostType ghost_type) {
AKANTU_DEBUG_IN();
UInt spatial_dimension = model->getSpatialDimension();
Mesh & mesh = model->getFEM().getMesh();
Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type);
Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type);
for(; it != last_type; ++it) {
switch(spatial_dimension) {
case 1: { assembleStiffnessMatrix<1>(current_position, *it, ghost_type); break; }
case 2: { assembleStiffnessMatrix<2>(current_position, *it, ghost_type); break; }
case 3: { assembleStiffnessMatrix<3>(current_position, *it, ghost_type); break; }
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
void Material::assembleStiffnessMatrix(Vector<Real> & current_position,
const ElementType & type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
SparseMatrix & K = const_cast<SparseMatrix &>(model->getStiffnessMatrix());
const Vector<Real> & shapes_derivatives = model->getFEM().getShapesDerivatives(type,ghost_type);
Vector<UInt> & elem_filter = element_filter(type, ghost_type);
Vector<Real> & strain_vect = strain(type, ghost_type);
UInt * elem_filter_val = elem_filter.storage();
UInt nb_element = elem_filter.getSize();
UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent();
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(type, ghost_type);
strain_vect.resize(nb_quadrature_points * nb_element);
model->getFEM().gradientOnQuadraturePoints(current_position, strain_vect,
dim, type, ghost_type, &elem_filter);
UInt tangent_size = getTangentStiffnessVoigtSize(dim);
Vector<Real> * tangent_stiffness_matrix =
new Vector<Real>(nb_element*nb_quadrature_points, tangent_size * tangent_size,
"tangent_stiffness_matrix");
computeTangentStiffness(type, *tangent_stiffness_matrix, ghost_type);
/// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$
UInt bt_d_b_size = dim * nb_nodes_per_element;
Vector<Real> * bt_d_b = new Vector<Real>(nb_element * nb_quadrature_points,
bt_d_b_size * bt_d_b_size,
"B^t*D*B");
UInt size_of_b = tangent_size * bt_d_b_size;
Real * B = new Real[size_of_b];
Real * Bt_D = new Real[size_of_b];
Real * Bt_D_B = bt_d_b->storage();
Real * D = tangent_stiffness_matrix->storage();
UInt offset_bt_d_b = bt_d_b_size * bt_d_b_size;
UInt offset_d = tangent_size * tangent_size;
for (UInt e = 0; e < nb_element; ++e) {
Real * shapes_derivatives_val =
shapes_derivatives.values + elem_filter_val[e]*size_of_shapes_derivatives*nb_quadrature_points;
for (UInt q = 0; q < nb_quadrature_points; ++q) {
transferBMatrixToSymVoigtBMatrix<dim>(shapes_derivatives_val, B, nb_nodes_per_element);
Math::matrixt_matrix(bt_d_b_size, tangent_size, tangent_size, B, D, Bt_D);
Math::matrix_matrix(bt_d_b_size, bt_d_b_size, tangent_size, Bt_D, B, Bt_D_B);
shapes_derivatives_val += size_of_shapes_derivatives;
D += offset_d;
Bt_D_B += offset_bt_d_b;
}
}
delete [] B;
delete [] Bt_D;
delete tangent_stiffness_matrix;
/// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$
Vector<Real> * K_e = new Vector<Real>(nb_element,
bt_d_b_size * bt_d_b_size,
"K_e");
model->getFEM().integrate(*bt_d_b, *K_e,
bt_d_b_size * bt_d_b_size,
type, ghost_type,
&elem_filter);
delete bt_d_b;
model->getFEM().assembleMatrix(*K_e, K, spatial_dimension, type, ghost_type, &elem_filter);
delete K_e;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void Material::computePotentialEnergy(ElementType el_type, GhostType ghost_type) {
AKANTU_DEBUG_IN();
if(ghost_type != _not_ghost) return;
Real * epot = potential_energy(el_type, ghost_type).storage();
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN;
computePotentialEnergy(strain_val, stress_val, epot);
epot++;
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void Material::computePotentialEnergyByElement() {
AKANTU_DEBUG_IN();
Mesh & mesh = model->getFEM().getMesh();
Mesh::type_iterator it = mesh.firstType(spatial_dimension);
Mesh::type_iterator last_type = mesh.lastType(spatial_dimension);
for(; it != last_type; ++it) {
if(!potential_energy.exists(*it, _not_ghost)) {
UInt nb_element = element_filter(*it, _not_ghost).getSize();
UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(*it, _not_ghost);
potential_energy.alloc(nb_element * nb_quadrature_points, 1,
*it, _not_ghost);
}
computePotentialEnergy(*it);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
Real Material::getPotentialEnergy() {
AKANTU_DEBUG_IN();
Real epot = 0.;
computePotentialEnergyByElement();
/// integrate the potential energy for each type of elements
Mesh & mesh = model->getFEM().getMesh();
Mesh::type_iterator it = mesh.firstType(spatial_dimension);
Mesh::type_iterator last_type = mesh.lastType(spatial_dimension);
for(; it != last_type; ++it) {
epot += model->getFEM().integrate(potential_energy(*it, _not_ghost), *it,
_not_ghost, &element_filter(*it, _not_ghost));
}
AKANTU_DEBUG_OUT();
return epot;
}
/* -------------------------------------------------------------------------- */
void Material::computeQuadraturePointsCoordinates(ByElementTypeReal & quadrature_points_coordinates) const {
AKANTU_DEBUG_IN();
const Mesh & mesh = model->getFEM().getMesh();
mesh.initByElementTypeVector(quadrature_points_coordinates, spatial_dimension, 0);
Vector<Real> nodes_coordinates(mesh.getNodes(), true);
nodes_coordinates += model->getDisplacement();
for(UInt gt = (UInt) _not_ghost; gt < (UInt) _casper; ++gt) {
GhostType ghost_type = (GhostType) gt;
Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type);
Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type);
for(; it != last_type; ++it) {
const Vector<UInt> & elem_filter = element_filter(*it, ghost_type);
UInt nb_element = elem_filter.getSize();
UInt nb_tot_quad = model->getFEM().getNbQuadraturePoints(*it, ghost_type) * nb_element;
Vector<Real> & quads = quadrature_points_coordinates(*it, ghost_type);
quads.resize(nb_tot_quad);
model->getFEM().interpolateOnQuadraturePoints(nodes_coordinates,
quads, spatial_dimension,
*it, ghost_type, &elem_filter);
}
}
AKANTU_DEBUG_OUT();
}
+/* -------------------------------------------------------------------------- */
+const Vector<Real> & Material::getVector(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) const {
+ std::stringstream sstr;
+ std::string ghost_id = "";
+ if (ghost_type == _ghost) ghost_id = ":ghost";
+ sstr << id << ":" << vect_id << ":" << type << ghost_id;
+
+ ID fvect_id = sstr.str();
+ try {
+ return Memory::getVector<Real>(fvect_id);
+ } catch(debug::Exception & e) {
+ AKANTU_EXCEPTION("The material " << name << "(" <<id << ") does not contain a vector " << vect_id << "(" << fvect_id << ") [" << e << "]");
+ }
+}
+
+/* -------------------------------------------------------------------------- */
+Real Material::getParam(const ID & param) const {
+ ID lparam = to_lower(param);
+ if(lparam == "rho") {
+ return rho;
+ }
+ AKANTU_EXCEPTION("No parameter named " << param << " in the material " << name << " (" << id << ")");
+}
+
+/* -------------------------------------------------------------------------- */
+void Material::setParam(const ID & param, Real value) {
+ ID lparam = to_lower(lparam);
+ if(lparam == "rho") {
+ rho = value;
+ }
+ AKANTU_EXCEPTION("No parameter named " << param << " in the material " << name << " (" << id << ")");
+}
+
+
/* -------------------------------------------------------------------------- */
void Material::printself(std::ostream & stream, int indent) const {
std::string space;
for(Int i = 0; i < indent; i++, space += AKANTU_INDENT);
stream << space << "Material [" << std::endl;
stream << space << " + id : " << id << std::endl;
stream << space << " + name : " << name << std::endl;
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
template void Material::initInternalVector<Real>(ByElementTypeVector<Real> & vect,
UInt nb_component,
ElementKind element_kind) const;
template void Material::initInternalVector<UInt>(ByElementTypeVector<UInt> & vect,
UInt nb_component,
ElementKind element_kind) const;
template void Material::initInternalVector<Int>(ByElementTypeVector<Int> & vect,
UInt nb_component,
ElementKind element_kind) const;
template void Material::resizeInternalVector<Real>(ByElementTypeVector<Real> & vect,
ElementKind element_kind) const;
template void Material::resizeInternalVector<UInt>(ByElementTypeVector<UInt> & vect,
ElementKind element_kind) const;
template void Material::resizeInternalVector<Int>(ByElementTypeVector<Int> & vect,
ElementKind element_kind) const;
__END_AKANTU__
diff --git a/src/model/solid_mechanics/material.hh b/src/model/solid_mechanics/material.hh
index 39711534a..9590859d5 100644
--- a/src/model/solid_mechanics/material.hh
+++ b/src/model/solid_mechanics/material.hh
@@ -1,432 +1,437 @@
/**
* @file material.hh
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Fri Jul 23 09:06:29 2010
*
* @brief Mother class for all materials
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "aka_memory.hh"
//#include "fem.hh"
//#include "mesh.hh"
#include "data_accessor.hh"
//#include "static_communicator.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MATERIAL_HH__
#define __AKANTU_MATERIAL_HH__
/* -------------------------------------------------------------------------- */
namespace akantu {
class Model;
class SolidMechanicsModel;
class CommunicationBuffer;
}
__BEGIN_AKANTU__
/**
* Interface of all materials
* Prerequisites for a new material
* - inherit from this class
* - implement the following methods:
* \code
* virtual Real getStableTimeStep(Real h, const Element & element = ElementNull);
*
* virtual void computeStress(ElementType el_type,
* GhostType ghost_type = _not_ghost);
*
* virtual void computeTangentStiffness(const ElementType & el_type,
* Vector<Real> & tangent_matrix,
* GhostType ghost_type = _not_ghost);
* \endcode
*
*/
class Material : protected Memory, public DataAccessor {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
Material(SolidMechanicsModel & model, const ID & id = "");
virtual ~Material();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// read properties
virtual bool setParam(const std::string & key, const std::string & value,
const ID & id);
/// initialize the material computed parameter
virtual void initMaterial();
/// compute the residual for this material
virtual void updateResidual(Vector<Real> & displacement,
GhostType ghost_type = _not_ghost);
void assembleResidual(GhostType ghost_type);
/// compute the residual for this material
virtual void computeStress(Vector<Real> & current_position,
GhostType ghost_type = _not_ghost);
/// compute the stiffness matrix
void assembleStiffnessMatrix(Vector<Real> & current_position,
GhostType ghost_type);
/// compute the stable time step for an element of size h
virtual Real getStableTimeStep(Real h, const Element & element = ElementNull) = 0;
/// compute the p-wave speed in the material
virtual Real getPushWaveSpeed() { AKANTU_DEBUG_TO_IMPLEMENT(); };
/// compute the s-wave speed in the material
virtual Real getShearWaveSpeed() { AKANTU_DEBUG_TO_IMPLEMENT(); };
/// add an element to the local mesh filter
inline UInt addElement(const ElementType & type,
UInt element,
const GhostType & ghost_type);
/// function to print the contain of the class
virtual void printself(std::ostream & stream, int indent = 0) const;
protected:
/// constitutive law
virtual void computeStress(ElementType el_type,
GhostType ghost_type = _not_ghost) = 0;
// /// constitutive law
// virtual void computeNonLocalStress(ElementType el_type,
// GhostType ghost_type = _not_ghost) {
// AKANTU_DEBUG_TO_IMPLEMENT();
// };
/// compute the tangent stiffness matrix
virtual void computeTangentStiffness(const ElementType & el_type,
Vector<Real> & tangent_matrix,
GhostType ghost_type = _not_ghost) = 0;
/// compute the potential energy
virtual void computePotentialEnergy(ElementType el_type,
GhostType ghost_type = _not_ghost);
template<UInt dim>
void assembleStiffnessMatrix(Vector<Real> & current_position,
const ElementType & type,
GhostType ghost_type);
/// transfer the B matrix to a Voigt notation B matrix
template<UInt dim>
inline void transferBMatrixToSymVoigtBMatrix(Real * B, Real * Bvoigt, UInt nb_nodes_per_element) const;
inline UInt getTangentStiffnessVoigtSize(UInt spatial_dimension) const;
/// compute the potential energy by element
void computePotentialEnergyByElement();
void computeQuadraturePointsCoordinates(ByElementTypeReal & quadrature_points_coordinates) const;
/* ------------------------------------------------------------------------ */
/* Function for all materials */
/* ------------------------------------------------------------------------ */
protected:
/// compute the potential energy for on element
inline void computePotentialEnergy(Real * F, Real * sigma, Real * epot);
public:
/// allocate an internal vector
template<typename T>
void initInternalVector(ByElementTypeVector<T> & vect,
UInt nb_component,
ElementKind element_kind = _ek_regular) const;
/// resize an internal vector
template<typename T>
void resizeInternalVector(ByElementTypeVector<T> & vect,
ElementKind element_kind = _ek_regular) const;
/* ------------------------------------------------------------------------ */
/* DataAccessor inherited members */
/* ------------------------------------------------------------------------ */
public:
virtual inline UInt getNbDataToPack(__attribute__((unused)) const Element & element,
__attribute__((unused)) SynchronizationTag tag) const {
return 0;
}
virtual inline UInt getNbDataToUnpack(__attribute__((unused)) const Element & element,
__attribute__((unused)) SynchronizationTag tag) const {
return 0;
}
virtual UInt getNbDataToPack(__attribute__((unused)) SynchronizationTag tag) const {
return 0;
}
virtual UInt getNbDataToUnpack(__attribute__((unused)) SynchronizationTag tag) const {
return 0;
}
virtual inline void packData(__attribute__((unused)) CommunicationBuffer & buffer,
__attribute__((unused)) const Element & element,
__attribute__((unused)) SynchronizationTag tag) const {
}
virtual void packData(__attribute__((unused)) CommunicationBuffer & buffer,
__attribute__((unused)) const UInt index,
__attribute__((unused)) SynchronizationTag tag) const {
}
virtual inline void unpackData(__attribute__((unused)) CommunicationBuffer & buffer,
__attribute__((unused)) const Element & element,
__attribute__((unused)) SynchronizationTag tag) {
}
virtual void unpackData(__attribute__((unused)) CommunicationBuffer & buffer,
__attribute__((unused)) const UInt index,
__attribute__((unused)) SynchronizationTag tag) {
}
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
AKANTU_GET_MACRO(Model, *model, const SolidMechanicsModel &)
AKANTU_GET_MACRO(ID, id, const ID &);
AKANTU_GET_MACRO(Rho, rho, Real);
AKANTU_SET_MACRO(Rho, rho, Real);
Real getPotentialEnergy();
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementFilter, element_filter, UInt);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Strain, strain, Real);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(PotentialEnergy, potential_energy, Real);
+ const Vector<Real> & getVector(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost) const;
+
+ virtual Real getParam(const ID & param) const;
+ virtual void setParam(const ID & param, Real value);
+
protected:
bool isInit() const { return is_init; }
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// id of the material
ID id;
/// material name
std::string name;
/// The model to witch the material belong
SolidMechanicsModel * model;
/// density : rho
Real rho;
/// stresses arrays ordered by element types
ByElementTypeReal stress;
/// strains arrays ordered by element types
ByElementTypeReal strain;
/// list of element handled by the material
ByElementTypeUInt element_filter;
/// is the vector for potential energy initialized
// bool potential_energy_vector;
/// potential energy by element
ByElementTypeReal potential_energy;
/// tell if using in non local mode or not
bool is_non_local;
/// spatial dimension
UInt spatial_dimension;
private:
/// boolean to know if the material has been initialized
bool is_init;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#if defined (AKANTU_INCLUDE_INLINE_IMPL)
# include "material_inline_impl.cc"
#endif
/// standard output stream operator
inline std::ostream & operator <<(std::ostream & stream, const Material & _this)
{
_this.printself(stream);
return stream;
}
__END_AKANTU__
/* -------------------------------------------------------------------------- */
/* Auto loop */
/* -------------------------------------------------------------------------- */
#define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN \
UInt nb_quadrature_points = this->model->getFEM().getNbQuadraturePoints(el_type, ghost_type); \
UInt size_strain = spatial_dimension * spatial_dimension; \
\
UInt nb_element = this->element_filter(el_type, ghost_type).getSize(); \
if (nb_element == 0) return; \
\
Vector<Real> & stress_tmp = this->stress(el_type, ghost_type); \
stress_tmp.resize(nb_element*nb_quadrature_points); \
Vector<Real> & strain_tmp = this->strain(el_type, ghost_type); \
\
Real * strain_val = strain_tmp.storage(); \
Real * stress_val = stress_tmp.storage(); \
\
for (UInt el = 0; el < nb_element; ++el) { \
for (UInt q = 0; q < nb_quadrature_points; ++q) { \
#define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END \
strain_val += size_strain; \
stress_val += size_strain; \
} \
} \
#define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent) \
UInt nb_quadrature_points = \
this->model->getFEM().getNbQuadraturePoints(el_type, ghost_type); \
UInt size_strain = spatial_dimension * spatial_dimension; \
\
UInt nb_element = this->element_filter(el_type, ghost_type).getSize(); \
if (nb_element == 0) return; \
\
Vector<Real> & strain_tmp = this->strain(el_type, ghost_type); \
\
Real * strain_val = strain_tmp.storage(); \
\
Real * tangent_val = tangent.values; \
UInt size_tangent = this->getTangentStiffnessVoigtSize(spatial_dimension); \
size_tangent *= size_tangent; \
\
for (UInt el = 0; el < nb_element; ++el) { \
for (UInt q = 0; q < nb_quadrature_points; ++q) { \
// Vector<Real> * stress_tmp = stress(el_type, ghost_type);
// stress_tmp->resize(nb_element*nb_quadrature_points);
// Real * stress_val = stress_tmp->values;
#define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END \
strain_val += size_strain; \
tangent_val += size_tangent; \
} \
}
/* -------------------------------------------------------------------------- */
/* Material list */
/* -------------------------------------------------------------------------- */
#define AKANTU_MATERIAL_LIST \
((elastic , MaterialElastic )) \
((elastic_orthotropic , MaterialElasticOrthotropic )) \
((viscoelastic , MaterialViscoElastic )) \
((elastic_caughey , MaterialElasticCaughey )) \
((neohookean , MaterialNeohookean )) \
((damage_linear , MaterialDamageLinear )) \
((marigo , MaterialMarigo )) \
((mazars , MaterialMazars )) \
((vreepeerlings , MaterialVreePeerlings )) \
((marigo_non_local , MaterialMarigoNonLocal )) \
((mazars_non_local , MaterialMazarsNonLocal )) \
((vreepeerlings_non_local, MaterialVreePeerlingsNonLocal)) \
((cohesive_bilinear , MaterialCohesiveBilinear )) \
((cohesive_linear , MaterialCohesiveLinear )) \
((cohesive_linear_extrinsic, MaterialCohesiveLinearExtrinsic )) \
((cohesive_linear_exponential_extrinsic, MaterialCohesiveLinearExponentialExtrinsic ))
#define INSTANSIATE_MATERIAL(mat_name) \
template class mat_name<1>; \
template class mat_name<2>; \
template class mat_name<3>
#if defined(__INTEL_COMPILER)
#pragma warning ( push )
/* warning #654: overloaded virtual function
"akantu::Material::computeStress" is only partially overridden in
class "akantu::Material*" */
#pragma warning ( disable : 654 )
#endif //defined(__INTEL_COMPILER)
/* -------------------------------------------------------------------------- */
// elastic materials
#include "material_elastic.hh"
#include "material_elastic_caughey.hh"
#include "material_viscoelastic.hh"
#include "material_neohookean.hh"
#include "material_elastic_orthotropic.hh"
// damage materials
#include "material_damage.hh"
#include "material_marigo.hh"
#include "material_mazars.hh"
#include "material_damage_linear.hh"
#include "material_vreepeerlings.hh"
#include "material_marigo_non_local.hh"
#include "material_mazars_non_local.hh"
#include "material_vreepeerlings_non_local.hh"
// cohesive materials
#include "material_cohesive.hh"
#include "material_cohesive_linear.hh"
#include "material_cohesive_bilinear.hh"
#include "material_cohesive_linear_extrinsic.hh"
#include "material_cohesive_linear_exponential_extrinsic.hh"
#if defined(__INTEL_COMPILER)
#pragma warning ( pop )
#endif //defined(__INTEL_COMPILER)
#endif /* __AKANTU_MATERIAL_HH__ */
diff --git a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.hh b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.hh
index 2791e6e6d..bfdb86df9 100644
--- a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.hh
+++ b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.hh
@@ -1,214 +1,214 @@
/**
* @file material_cohesive.hh
* @author Marco Vocialta <marco.vocialta@epfl.ch>
* @date Tue Feb 7 17:50:23 2012
*
* @brief Specialization of the material class for cohesive elements
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "material.hh"
#include "aka_common.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MATERIAL_COHESIVE_HH__
#define __AKANTU_MATERIAL_COHESIVE_HH__
/* -------------------------------------------------------------------------- */
namespace akantu {
class SolidMechanicsModelCohesive;
}
__BEGIN_AKANTU__
class MaterialCohesive : public Material {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
typedef FEMTemplate< IntegratorCohesive<IntegratorGauss>,
ShapeCohesive<ShapeLagrange> > MyFEMCohesiveType;
public:
MaterialCohesive(SolidMechanicsModel& model, const ID & id = "");
virtual ~MaterialCohesive();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// read properties
virtual bool setParam(const std::string & key, const std::string & value,
const ID & id);
/// initialize the material computed parameter
virtual void initMaterial();
/// resize vectors for new cohesive elements
virtual void resizeCohesiveVectors();
/// compute the residual for this material
virtual void updateResidual(Vector<Real> & current_position,
GhostType ghost_type = _not_ghost);
/// compute the stable time step for an element of size h
virtual Real getStableTimeStep(__attribute__((unused)) Real h,
__attribute__((unused)) const Element & element = ElementNull) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
// /// add an element to the local mesh filter
// __aka_inline__ void addElement(const ElementType & type,
// UInt element,
// const GhostType & ghost_type);
/// function to print the contain of the class
virtual void printself(std::ostream & stream, int indent = 0) const;
/// check stress for cohesive elements' insertion
virtual void checkInsertion(Vector<UInt> & facet_insertion);
protected:
/// constitutive law
virtual void computeStress(__attribute__((unused)) ElementType el_type,
__attribute__((unused)) GhostType ghost_type = _not_ghost) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/// compute the tangent stiffness matrix
virtual void computeTangentStiffness(__attribute__((unused)) const ElementType & el_type,
__attribute__((unused)) Vector<Real> & tangent_matrix,
__attribute__((unused)) GhostType ghost_type = _not_ghost) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
void computeNormal(const Vector<Real> & position,
Vector<Real> & normal,
ElementType type,
GhostType ghost_type);
void computeOpening(const Vector<Real> & displacement,
Vector<Real> & normal,
ElementType type,
GhostType ghost_type);
template<ElementType type>
void computeNormal(const Vector<Real> & position,
Vector<Real> & normal,
GhostType ghost_type);
/// assemble residual
void assembleResidual(GhostType ghost_type = _not_ghost);
/// compute tractions (including normals and openings)
void computeTraction(GhostType ghost_type = _not_ghost);
/// constitutive law
virtual void computeTraction(const Vector<Real> & normal,
ElementType el_type,
GhostType ghost_type = _not_ghost) = 0;
/// compute reversible and total energies by element
void computeEnergies();
/// compute effective stress norm for insertion check
- virtual Real computeEffectiveNorm(const types::Matrix & stress,
- const types::RVector & normal,
- const types::RVector & tangent){
+ virtual Real computeEffectiveNorm(__attribute__((unused)) const types::Matrix & stress,
+ __attribute__((unused)) const types::RVector & normal,
+ __attribute__((unused)) const types::RVector & tangent){
AKANTU_DEBUG_TO_IMPLEMENT();
};
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// compute reversible energy
Real getReversibleEnergy();
/// compute dissipated energy
Real getDissipatedEnergy();
/// get sigma_c
AKANTU_GET_MACRO(SigmaC, sigma_c, Real);
/// get rand
AKANTU_GET_MACRO(RandFactor, rand, Real);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
/// reversible energy by quadrature point
ByElementTypeReal reversible_energy;
/// total energy by quadrature point
ByElementTypeReal total_energy;
/// traction in all elements and quadrature points (previous time step)
ByElementTypeReal tractions_old;
/// opening in all elements and quadrature points (previous time step)
ByElementTypeReal opening_old;
protected:
/// traction in all elements and quadrature points
ByElementTypeReal tractions;
/// opening in all elements and quadrature points
ByElementTypeReal opening;
/// Link to the cohesive fem object in the model
MyFEMCohesiveType * fem_cohesive;
/// critical stress
Real sigma_c;
/// randomness factor
Real rand;
/// vector to store stresses on facets for element insertions
Vector<Real> sigma_insertion;
/// pointer to the solid mechanics model for cohesive elements
SolidMechanicsModelCohesive * model;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#include "material_cohesive_inline_impl.cc"
__END_AKANTU__
#endif /* __AKANTU_MATERIAL_COHESIVE_HH__ */
diff --git a/src/model/solid_mechanics/materials/material_elastic.cc b/src/model/solid_mechanics/materials/material_elastic.cc
index 9703abc2e..f2e4e82ca 100644
--- a/src/model/solid_mechanics/materials/material_elastic.cc
+++ b/src/model/solid_mechanics/materials/material_elastic.cc
@@ -1,166 +1,191 @@
/**
* @file material_elastic.cc
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Tue Jul 27 11:53:52 2010
*
* @brief Specialization of the material class for the elastic material
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "material_elastic.hh"
#include "solid_mechanics_model.hh"
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
MaterialElastic<spatial_dimension>::MaterialElastic(SolidMechanicsModel & model, const ID & id) :
Material(model, id) {
AKANTU_DEBUG_IN();
E = 0;
nu = 1./2.;
plane_stress = false;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialElastic<spatial_dimension>::initMaterial() {
AKANTU_DEBUG_IN();
Material::initMaterial();
lambda = nu * E / ((1 + nu) * (1 - 2*nu));
mu = E / (2 * (1 + nu));
if(plane_stress) {
//lambda = 2 * lambda * mu / (lambda + 2 * mu);
lambda = nu * E / ((1 + nu)*(1 - nu));
}
kpa = lambda + 2./3. * mu;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialElastic<spatial_dimension>::computeStress(ElementType el_type, GhostType ghost_type) {
AKANTU_DEBUG_IN();
Real F[3*3];
Real sigma[3*3];
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN;
memset(F, 0, 3 * 3 * sizeof(Real));
for (UInt i = 0; i < spatial_dimension; ++i)
for (UInt j = 0; j < spatial_dimension; ++j)
F[3*i + j] = strain_val[spatial_dimension * i + j];
// for (UInt i = 0; i < spatial_dimension; ++i) F[i*3 + i] -= 1;
computeStress(F, sigma);
for (UInt i = 0; i < spatial_dimension; ++i)
for (UInt j = 0; j < spatial_dimension; ++j)
stress_val[spatial_dimension*i + j] = sigma[3 * i + j];
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
-void MaterialElastic<spatial_dimension>::computeTangentStiffness(const ElementType & el_type,
+void MaterialElastic<spatial_dimension>::computeTangentStiffness(__attribute__((unused)) const ElementType & el_type,
Vector<Real> & tangent_matrix,
- GhostType ghost_type) {
+ __attribute__((unused)) GhostType ghost_type) {
AKANTU_DEBUG_IN();
Real * tangent_val = tangent_matrix.values;
UInt offset_tangent = tangent_matrix.getNbComponent();
UInt nb_quads = tangent_matrix.getSize();
if (nb_quads == 0) return;
memset(tangent_val, 0, offset_tangent * nb_quads * sizeof(Real));
for (UInt q = 0; q < nb_quads; ++q, tangent_val += offset_tangent) {
computeTangentStiffness(tangent_val);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
bool MaterialElastic<spatial_dimension>::setParam(const std::string & key, const std::string & value,
const ID & id) {
std::stringstream sstr(value);
if(key == "E") { sstr >> E; }
else if(key == "nu") { sstr >> nu; }
else if(key == "Plane_Stress") { sstr >> plane_stress; }
else { return Material::setParam(key, value, id); }
return true;
}
+/* -------------------------------------------------------------------------- */
+template<UInt spatial_dimension>
+Real MaterialElastic<spatial_dimension>::getParam(const ID & param) const {
+ ID key = to_lower(param);
+ if(key == "e") { return E; }
+ else if(key == "nu") { return nu; }
+ else if(key == "lambda") { return nu; }
+ else if(key == "mu") { return mu; }
+ else if(key == "kapa") { return kpa; }
+ else return Material::getParam(param);
+}
+
+/* -------------------------------------------------------------------------- */
+template<UInt spatial_dimension>
+void MaterialElastic<spatial_dimension>::setParam(const ID & param, Real value) {
+ ID key = to_lower(param);
+ if(key == "e") { E = value; }
+ else if(key == "nu") { nu = value; }
+ else if(key == "lambda") { nu = value; }
+ else if(key == "mu") { mu = value; }
+ else if(key == "kapa") { kpa = value; }
+ else Material::setParam(param, value);
+}
+
+
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
Real MaterialElastic<spatial_dimension>::getPushWaveSpeed() {
return sqrt((lambda + 2*mu)/rho);
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
Real MaterialElastic<spatial_dimension>::getShearWaveSpeed() {
return sqrt(mu/rho);
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialElastic<spatial_dimension>::printself(std::ostream & stream,
int indent) const {
std::string space;
for(Int i = 0; i < indent; i++, space += AKANTU_INDENT);
stream << space << "Material<_elastic> [" << std::endl;
if(!plane_stress)
stream << space << " + Plane strain" << std::endl;
else
stream << space << " + Plane stress" << std::endl;
stream << space << " + density : " << rho << std::endl;
stream << space << " + Young's modulus : " << E << std::endl;
stream << space << " + Poisson's ratio : " << nu << std::endl;
if(this->isInit()) {
stream << space << " + First Lamé coefficient : " << lambda << std::endl;
stream << space << " + Second Lamé coefficient : " << mu << std::endl;
stream << space << " + Bulk coefficient : " << kpa << std::endl;
}
Material::printself(stream, indent + 1);
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
INSTANSIATE_MATERIAL(MaterialElastic);
__END_AKANTU__
diff --git a/src/model/solid_mechanics/materials/material_elastic.hh b/src/model/solid_mechanics/materials/material_elastic.hh
index 9fb179e6d..b6c5b4d72 100644
--- a/src/model/solid_mechanics/materials/material_elastic.hh
+++ b/src/model/solid_mechanics/materials/material_elastic.hh
@@ -1,147 +1,151 @@
/**
* @file material_elastic.hh
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Thu Jul 29 15:00:59 2010
*
* @brief Material isotropic elastic
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "material.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MATERIAL_ELASTIC_HH__
#define __AKANTU_MATERIAL_ELASTIC_HH__
__BEGIN_AKANTU__
/**
* Material elastic isotropic
*
* parameters in the material files :
* - rho : density (default: 0)
* - E : Young's modulus (default: 0)
* - nu : Poisson's ratio (default: 1/2)
* - Plane_Stress : if 0: plane strain, else: plane stress (default: 0)
*/
template<UInt spatial_dimension>
class MaterialElastic : public virtual Material {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
MaterialElastic(SolidMechanicsModel & model, const ID & id = "");
virtual ~MaterialElastic() {};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
virtual void initMaterial();
virtual bool setParam(const std::string & key, const std::string & value,
const ID & id);
/// constitutive law for all element of a type
virtual void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost);
/// compute the tangent stiffness matrix for an element type
void computeTangentStiffness(const ElementType & el_type,
Vector<Real> & tangent_matrix,
GhostType ghost_type = _not_ghost);
/// compute the p-wave speed in the material
virtual Real getPushWaveSpeed();
/// compute the s-wave speed in the material
virtual Real getShearWaveSpeed();
/// function to print the containt of the class
virtual void printself(std::ostream & stream, int indent = 0) const;
protected:
/// constitutive law for a given quadrature point
inline void computeStress(Real * F, Real * sigma);
inline void computeStress(types::Matrix & grad_u,
types::Matrix & sigma);
// /// compute the tangent stiffness matrix for an element
void computeTangentStiffness(Real * tangent);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// get the stable time step
inline Real getStableTimeStep(Real h, const Element & element);
AKANTU_GET_MACRO(E, E, Real);
AKANTU_GET_MACRO(Nu, nu, Real);
AKANTU_GET_MACRO(Mu, mu, Real);
AKANTU_GET_MACRO(Lambda, lambda, Real);
AKANTU_GET_MACRO(Kpa, kpa, Real);
AKANTU_GET_MACRO(PlaneStress, plane_stress, bool);
AKANTU_SET_MACRO(E, E, Real);
AKANTU_SET_MACRO(Nu, nu, Real);
+
+ virtual Real getParam(const ID & param) const;
+ virtual void setParam(const ID & param, Real value);
+
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// the young modulus
Real E;
/// Poisson coefficient
Real nu;
/// First Lamé coefficient
Real lambda;
/// Second Lamé coefficient (shear modulus)
Real mu;
/// Bulk modulus
Real kpa;
/// Plane stress or plane strain
bool plane_stress;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#include "material_elastic_inline_impl.cc"
__END_AKANTU__
#endif /* __AKANTU_MATERIAL_ELASTIC_HH__ */
diff --git a/src/model/solid_mechanics/materials/material_elastic_caughey.cc b/src/model/solid_mechanics/materials/material_elastic_caughey.cc
index fe1e613f5..e5c148814 100644
--- a/src/model/solid_mechanics/materials/material_elastic_caughey.cc
+++ b/src/model/solid_mechanics/materials/material_elastic_caughey.cc
@@ -1,174 +1,191 @@
/**
* @file material_elastic_caughey.cc
* @author David Kammer <david.kammer@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Wed May 4 15:25:52 2011
*
* @brief Special. of the material class for the caughey viscoelastic material
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "material_elastic_caughey.hh"
#include "solid_mechanics_model.hh"
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
MaterialElasticCaughey<spatial_dimension>::MaterialElasticCaughey(SolidMechanicsModel & model,
const ID & id) :
Material(model, id), MaterialElastic<spatial_dimension>(model, id),
stress_viscosity("stress_viscosity", id),
stress_elastic("stress_elastic", id) {
AKANTU_DEBUG_IN();
alpha = 0.;
this->initInternalVector(this->stress_viscosity, spatial_dimension * spatial_dimension);
this->initInternalVector(this->stress_elastic , spatial_dimension * spatial_dimension);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialElasticCaughey<spatial_dimension>::initMaterial() {
AKANTU_DEBUG_IN();
MaterialElastic<spatial_dimension>::initMaterial();
this->resizeInternalVector(this->stress_viscosity);
this->resizeInternalVector(this->stress_elastic );
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialElasticCaughey<spatial_dimension>::computeStress(ElementType el_type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
Vector<UInt> & elem_filter = this->element_filter (el_type, ghost_type);
Vector<Real> & stress_visc = stress_viscosity(el_type, ghost_type);
Vector<Real> & stress_el = stress_elastic (el_type, ghost_type);
MaterialElastic<spatial_dimension>::computeStress(el_type, ghost_type);
Vector<Real> & velocity = this->model->getVelocity();
UInt nb_elem = this->element_filter(el_type, ghost_type).getSize();
UInt nb_quad_points_per_elem =
this->model->getFEM().getNbQuadraturePoints(el_type, ghost_type);
UInt nb_quad_points = nb_quad_points_per_elem * nb_elem;
Vector<Real> strain_rate(nb_quad_points,
spatial_dimension * spatial_dimension,
"strain_rate");
this->model->getFEM().gradientOnQuadraturePoints(velocity, strain_rate,
spatial_dimension,
el_type, ghost_type, &elem_filter);
Real F[3*3];
Real sigma[3*3];
Real * strain_rate_val = strain_rate.storage();
Real * stress_visc_val = stress_visc.storage();
Real * stress_el_val = stress_el.storage();
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN;
memset(F, 0, 3 * 3 * sizeof(Real));
for (UInt i = 0; i < spatial_dimension; ++i)
for (UInt j = 0; j < spatial_dimension; ++j)
F[3*i + j] = strain_rate_val[spatial_dimension * i + j];
MaterialElastic<spatial_dimension>::computeStress(F, sigma);
for (UInt i = 0; i < spatial_dimension; ++i)
for (UInt j = 0; j < spatial_dimension; ++j) {
stress_visc_val[spatial_dimension * i + j] = alpha * sigma[3 * i + j];
stress_el_val [spatial_dimension * i + j] = stress_val[spatial_dimension*i + j];
stress_val [spatial_dimension * i + j] += stress_visc_val[spatial_dimension*i + j];
}
strain_rate_val += spatial_dimension * spatial_dimension;
stress_visc_val += spatial_dimension * spatial_dimension;
stress_el_val += spatial_dimension * spatial_dimension;
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialElasticCaughey<spatial_dimension>::computePotentialEnergy(ElementType el_type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
if(ghost_type != _not_ghost) return;
Vector<Real> & stress_el = stress_elastic(el_type, ghost_type);
Real * stress_el_val = stress_el.storage();
Real * epot = this->potential_energy(el_type, ghost_type).storage();
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN;
MaterialElastic<spatial_dimension>::computePotentialEnergy(strain_val, stress_el_val, epot);
epot++;
stress_el_val += spatial_dimension*spatial_dimension;
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
bool MaterialElasticCaughey<spatial_dimension>::setParam(const std::string & key,
const std::string & value,
const ID & id) {
std::stringstream sstr(value);
if(key == "alpha") { sstr >> alpha; }
else { return MaterialElastic<spatial_dimension>::setParam(key, value, id); }
return true;
}
+/* -------------------------------------------------------------------------- */
+template<UInt spatial_dimension>
+Real MaterialElasticCaughey<spatial_dimension>::getParam(const ID & param) const {
+ ID key = to_lower(param);
+ if(key == "alpha") { return alpha; }
+ else return Material::getParam(param);
+}
+
+/* -------------------------------------------------------------------------- */
+template<UInt spatial_dimension>
+void MaterialElasticCaughey<spatial_dimension>::setParam(const ID & param, Real value) {
+ ID key = to_lower(param);
+ if(key == "alpha") { alpha = value; }
+ else Material::setParam(param, value);
+}
+
+
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialElasticCaughey<spatial_dimension>::printself(std::ostream & stream,
int indent) const {
std::string space;
for(Int i = 0; i < indent; i++, space += AKANTU_INDENT);
stream << space << "MaterialElasticCaughey [" << std::endl;
MaterialElastic<spatial_dimension>::printself(stream, indent + 1);
stream << space << " + artifical viscous ratio : " << alpha << std::endl;
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
INSTANSIATE_MATERIAL(MaterialElasticCaughey);
__END_AKANTU__
diff --git a/src/model/solid_mechanics/materials/material_elastic_caughey.hh b/src/model/solid_mechanics/materials/material_elastic_caughey.hh
index 2ed74fc2e..5675661fd 100644
--- a/src/model/solid_mechanics/materials/material_elastic_caughey.hh
+++ b/src/model/solid_mechanics/materials/material_elastic_caughey.hh
@@ -1,118 +1,122 @@
/**
* @file material_elastic_caughey.hh
* @author David Kammer <david.kammer@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Wed May 4 15:16:59 2011
*
* @brief Material isotropic viscoelastic (according to the Caughey condition)
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "material.hh"
#include "material_elastic.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MATERIAL_ELASTIC_CAUGHEY_HH__
#define __AKANTU_MATERIAL_ELASTIC_CAUGHEY_HH__
__BEGIN_AKANTU__
/**
* Material viscoelastic (caughey condition) isotropic
*
* parameters in the material files :
* - rho : density (default: 0)
* - E : Young's modulus (default: 0)
* - nu : Poisson's ratio (default: 1/2)
* - Plane_Stress : if 0: plane strain, else: plane stress (default: 0)
* - alpha : viscous ratio
*/
template<UInt spatial_dimension>
class MaterialElasticCaughey : public MaterialElastic<spatial_dimension> {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
MaterialElasticCaughey(SolidMechanicsModel & model, const ID & id = "");
virtual ~MaterialElasticCaughey() {};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
void initMaterial();
virtual bool setParam(const std::string & key, const std::string & value,
const ID & id);
/// constitutive law for all element of a type
void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost);
/// compute the potential energy for all elements
virtual void computePotentialEnergy(ElementType el_type, GhostType ghost_type = _not_ghost);
/// function to print the containt of the class
virtual void printself(std::ostream & stream, int indent = 0) const;
protected:
/// constitutive law for a given quadrature point
//inline void computeStress(Real * F, Real * sigma);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(StressViscosity, stress_viscosity, Real);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(StressElastic, stress_elastic, Real);
AKANTU_GET_MACRO(Alpha, alpha, Real);
AKANTU_SET_MACRO(Alpha, alpha, Real);
+ virtual Real getParam(const ID & param) const;
+ virtual void setParam(const ID & param, Real value);
+
+
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
/// stress due to viscosity
ByElementTypeReal stress_viscosity;
/// stress due to elasticity
ByElementTypeReal stress_elastic;
/// viscous ratio
Real alpha;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
//#include "material_elastic_caughey_inline_impl.cc"
__END_AKANTU__
#endif /* __AKANTU_MATERIAL_ELASTIC_CAUGHEY_HH__ */
diff --git a/src/model/solid_mechanics/materials/material_viscoelastic.cc b/src/model/solid_mechanics/materials/material_viscoelastic.cc
index 6890ff6ca..939ef007e 100644
--- a/src/model/solid_mechanics/materials/material_viscoelastic.cc
+++ b/src/model/solid_mechanics/materials/material_viscoelastic.cc
@@ -1,174 +1,176 @@
/**
* @file material_viscoelastic.hh
* @author Vlad Yastrebov <vladislav.yastrebov@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author David Kammer <david.kammer@epfl.ch>
* @date Thu Feb 7 2012
*
* @brief Material Visco-elastic
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "material_viscoelastic.hh"
#include "solid_mechanics_model.hh"
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
MaterialViscoElastic<spatial_dimension>::MaterialViscoElastic(SolidMechanicsModel & model, const ID & id) :
Material(model, id), MaterialElastic<spatial_dimension>(model, id),
stress_dev("stress_dev", id),
history_integral("history_integral", id) {
AKANTU_DEBUG_IN();
eta = 1.;
E_inf = 1.;
Ev = 1.;
UInt stress_size = spatial_dimension * spatial_dimension;
this->initInternalVector(this->stress_dev, stress_size);
this->initInternalVector(this->history_integral, stress_size);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialViscoElastic<spatial_dimension>::initMaterial() {
AKANTU_DEBUG_IN();
E_inf = this->E - this->Ev;
MaterialElastic<spatial_dimension>::initMaterial();
this->resizeInternalVector(this->stress_dev);
this->resizeInternalVector(this->history_integral);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialViscoElastic<spatial_dimension>::computeStress(ElementType el_type, GhostType ghost_type) {
AKANTU_DEBUG_IN();
- Real tau = eta / Ev;
+ Real tau = 0.;
+ // if(std::abs(Ev) > std::numeric_limits<Real>::epsilon())
+ tau = eta / Ev;
Vector<Real> & stress_dev_vect = stress_dev(el_type, ghost_type);
Vector<Real> & history_int_vect = history_integral(el_type, ghost_type);
Vector<Real>::iterator<types::Matrix> stress_d = stress_dev_vect.begin(spatial_dimension, spatial_dimension);
Vector<Real>::iterator<types::Matrix> history_int = history_int_vect.begin(spatial_dimension, spatial_dimension);
types::Matrix e(spatial_dimension, spatial_dimension);
types::Matrix s(spatial_dimension, spatial_dimension);
types::Matrix theta_sp(spatial_dimension, spatial_dimension);
Real dt = this->model->getTimeStep();
Real exp_dt_tau = exp( -dt/tau );
Real exp_dt_tau_2 = exp( -.5*dt/tau );
Real plane_stress_coeff = this->nu / (this->nu - 1);
/// Loop on all quadrature points
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN;
types::Matrix sigma(stress_val, spatial_dimension, spatial_dimension);
types::Matrix grad_u(strain_val, spatial_dimension, spatial_dimension);
types::Matrix & dev_s = *stress_d;
types::Matrix & h = *history_int;
e.clear();
s.clear();
sigma.clear();
/// Compute the first invariant of strain
Real Theta = grad_u.trace();
///\todo correct the trace in case of plane stress
if(spatial_dimension == 2 && this->plane_stress == true)
Theta += plane_stress_coeff * (grad_u(0,0) + grad_u(1,1));
for (UInt i = 0; i < spatial_dimension; ++i)
for (UInt j = 0; j < spatial_dimension; ++j) {
s(i, j) = 2 * this->mu * (.5 * (grad_u(i,j) + grad_u(j,i)) - 1./3. * Theta *(i == j));
}
/// Update the internal variable
for (UInt i = 0; i < spatial_dimension; ++i)
for (UInt j = 0; j < spatial_dimension; ++j) {
h(i, j) = exp_dt_tau * h(i, j) + exp_dt_tau_2 * (s(i, j) - dev_s(i, j));
dev_s(i, j) = s(i, j);
}
types::Matrix U_rond_prim(spatial_dimension, spatial_dimension);
U_rond_prim.eye(this->kpa * Theta);
Real gamma_inf = E_inf / this->E;
Real gamma_v = Ev / this->E;
for (UInt i = 0; i < spatial_dimension; ++i)
for (UInt j = 0; j < spatial_dimension; ++j)
sigma(i, j) = U_rond_prim(i,j) + gamma_inf * s(i, j) + gamma_v * h(i, j);
/// Save the deviator of stress
++stress_d;
++history_int;
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
bool MaterialViscoElastic<spatial_dimension>::setParam(const std::string & key, const std::string & value,
const ID & id) {
std::stringstream sstr(value);
if(key == "eta") { sstr >> eta; }
else if(key == "Ev") { sstr >> Ev; }
else { return MaterialElastic<spatial_dimension>::setParam(key, value, id); }
return true;
}
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
void MaterialViscoElastic<spatial_dimension>::printself(std::ostream & stream, int indent) const {
std::string space;
for(Int i = 0; i < indent; i++, space += AKANTU_INDENT);
stream << space << "MaterialViscoElastic [" << std::endl;
MaterialElastic<spatial_dimension>::printself(stream, indent + 1);
stream << space << " + Eta : " << eta << std::endl;
stream << space << " + Ev : " << Ev << std::endl;
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
INSTANSIATE_MATERIAL(MaterialViscoElastic);
__END_AKANTU__
diff --git a/src/model/solid_mechanics/materials/material_viscoelastic.hh b/src/model/solid_mechanics/materials/material_viscoelastic.hh
index 3d9e8d792..61ff5cb92 100644
--- a/src/model/solid_mechanics/materials/material_viscoelastic.hh
+++ b/src/model/solid_mechanics/materials/material_viscoelastic.hh
@@ -1,139 +1,139 @@
/**
* @file material_viscoelastic.hh
* @author Vlad Yastrebov <vladislav.yastrebov@epfl.ch>
* @date Thu Feb 7 2012
*
* @brief Material Visco-elastic, based on Standard Solid rheological model, see
* [] J.C. Simo, T.J.R. Hughes, "Computational Inelasticity", Springer (1998),
* see Sections 10.2 and 10.3
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "material.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MATERIAL_VISCOELASTIC_HH__
#define __AKANTU_MATERIAL_VISCOELASTIC_HH__
__BEGIN_AKANTU__
/**
* Material viscoelastic
*
*
* @verbatim
E_\inf
------|\/\/\|------
| |
---| |---
| |
----|\/\/\|--[|----
E_v \eta
@endverbatim
*
* parameters in the material files :
* - E : Initial Young's modulus @f[ E = E_\inf + E_v @f]
* - eta : viscosity
* - Ev : stiffness of the viscous element
*/
template<UInt spatial_dimension>
class MaterialViscoElastic : public MaterialElastic<spatial_dimension> {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
MaterialViscoElastic(SolidMechanicsModel & model, const ID & id = "");
virtual ~MaterialViscoElastic() {};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
void initMaterial();
virtual bool setParam(const std::string & key, const std::string & value,
const ID & id);
/// constitutive law for all element of a type
void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost);
/// function to print the containt of the class
virtual void printself(std::ostream & stream, int indent = 0) const;
protected:
/// constitutive law for a given quadrature point
//__aka_inline__ void computeStress(Real * F, Real * sigma);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(HistoryIntegral, history_integral, Real);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(DeviatoricStress, stress_dev, Real);
- AKANTU_GET_MACRO(Ev, Ev, const Real);
+ AKANTU_GET_MACRO(Ev, Ev, Real);
AKANTU_SET_MACRO(Ev, Ev, Real);
- AKANTU_GET_MACRO(Einf, E_inf, const Real);
+ AKANTU_GET_MACRO(Einf, E_inf, Real);
// AKANTU_SET_MACRO(Einf, E_inf, Real); // Einf is computed in initMaterial based on E and Ev
- AKANTU_GET_MACRO(Eta, eta, const Real);
+ AKANTU_GET_MACRO(Eta, eta, Real);
AKANTU_SET_MACRO(Eta, eta, Real);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
/// viscosity, viscous elastic modulus
Real eta, Ev, E_inf;
/// history of deviatoric stress
ByElementTypeReal stress_dev;
/// Internal variable: history integral
ByElementTypeReal history_integral;
};
/* -------------------------------------------------------------------------- */
/* __aka_inline__ functions */
/* -------------------------------------------------------------------------- */
//#include "material_elastic_caughey_inline_impl.cc"
/* -------------------------------------------------------------------------- */
/// standard output stream operator
/*
inline std::ostream & operator <<(std::ostream & stream, const MaterialViscoElastic & _this)
{
_this.printself(stream);
return stream;
}
*/
__END_AKANTU__
#endif /* __AKANTU_MATERIAL_VISCOELASTIC_HH__ */
diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc
index 41d7d97d6..8d76db534 100644
--- a/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc
+++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc
@@ -1,675 +1,673 @@
/**
* @file solid_mechanics_model_cohesive.cc
* @author Marco Vocialta <marco.vocialta@epfl.ch>
* @date Thu Apr 19 10:42:11 2012
*
* @brief Solid mechanics model for cohesive elements
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include <cstdlib>
#include <algorithm>
#include <ctime>
#include "solid_mechanics_model_cohesive.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
SolidMechanicsModelCohesive::SolidMechanicsModelCohesive(Mesh & mesh,
UInt dim,
const ID & id,
const MemoryID & memory_id)
: SolidMechanicsModel(mesh, dim, id, memory_id),
mesh_facets(mesh.getSpatialDimension(),
mesh.getNodes().getID(),
id, memory_id) {
AKANTU_DEBUG_IN();
registerFEMObject<MyFEMCohesiveType>("CohesiveFEM", mesh, spatial_dimension);
MeshUtils::buildAllFacets(mesh, mesh_facets);
/// assign cohesive type
Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1);
Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1);
for (; it != last; ++it) {
const Vector<UInt> & connectivity = mesh_facets.getConnectivity(*it);
if (connectivity.getSize() != 0) {
type_facet = *it;
break;
}
}
if (type_facet == _segment_2) type_cohesive = _cohesive_2d_4;
else if (type_facet == _segment_3) type_cohesive = _cohesive_2d_6;
mesh.addConnectivityType(type_cohesive);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::initFull(std::string material_file,
AnalysisMethod method) {
SolidMechanicsModel::initFull(material_file, method);
if(material_file != "") {
initCohesiveMaterial();
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::initCohesiveMaterial() {
AKANTU_DEBUG_IN();
/// find the cohesive index
cohesive_index = 0;
- MaterialCohesive * mat_cohesive;
- while ((mat_cohesive =
- dynamic_cast<MaterialCohesive *>(materials[cohesive_index])) == NULL
+ while ((dynamic_cast<MaterialCohesive *>(materials[cohesive_index]) == NULL)
&& cohesive_index <= materials.size())
++cohesive_index;
AKANTU_DEBUG_ASSERT(cohesive_index != materials.size(),
"No cohesive materials in the material input file");
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Initialize the model,basically it pre-compute the shapes, shapes derivatives
* and jacobian
*
*/
void SolidMechanicsModelCohesive::initModel() {
SolidMechanicsModel::initModel();
getFEM("CohesiveFEM").initShapeFunctions(_not_ghost);
getFEM("CohesiveFEM").initShapeFunctions(_ghost);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::initExtrinsic() {
AKANTU_DEBUG_IN();
MaterialCohesive * mat_cohesive
= dynamic_cast<MaterialCohesive*>(materials[cohesive_index]);
UInt nb_facet = mesh_facets.getNbElement(type_facet);
sigma_lim.resize(nb_facet);
const Real sigma_c = mat_cohesive->getSigmaC();
const Real rand = mat_cohesive->getRandFactor();
std::srand(time(NULL));
const Vector<Vector<Element> > & element_to_facet
= mesh_facets.getElementToSubelement(type_facet);
facets_check.resize(nb_facet);
for (UInt f = 0; f < nb_facet; ++f) {
if (element_to_facet(f)(1) != ElementNull) {
facets_check(f) = true;
sigma_lim(f) = sigma_c * (1 + std::rand()/(Real)RAND_MAX * rand);
}
else {
facets_check(f) = false;
sigma_lim(f) = std::numeric_limits<Real>::max();
}
}
registerFEMObject<MyFEMType>("FacetsFEM", mesh_facets, spatial_dimension-1);
getFEM("FacetsFEM").initShapeFunctions();
getFEM("FacetsFEM").computeNormalsOnControlPoints();
/// THIS HAS TO BE CHANGED:
const Vector<Real> & normals = getFEM("FacetsFEM").getNormalsOnQuadPoints(type_facet);
Vector<Real>::const_iterator<types::RVector> normal_it =
normals.begin(spatial_dimension);
tangents.resize(normals.getSize());
tangents.extendComponents(spatial_dimension);
Vector<Real>::iterator<types::RVector> tangent_it =
tangents.begin(spatial_dimension);
for (UInt i = 0; i < normals.getSize(); ++i, ++normal_it, ++tangent_it) {
Math::normal2( (*normal_it).storage(), (*tangent_it).storage() );
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::checkCohesiveStress() {
AKANTU_DEBUG_IN();
Vector<UInt> facet_insertion;
MaterialCohesive * mat_cohesive
= dynamic_cast<MaterialCohesive*>(materials[cohesive_index]);
mat_cohesive->checkInsertion(facet_insertion);
if (facet_insertion.getSize() != 0)
insertCohesiveElements(facet_insertion);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::insertCohesiveElements(const Vector<UInt> & facet_insertion) {
AKANTU_DEBUG_IN();
if(facet_insertion.getSize() == 0) return;
Vector<UInt> doubled_nodes(0, 2);
Vector<UInt> doubled_facets(0, 2);
/// update mesh
for (UInt f = 0; f < facet_insertion.getSize(); ++f) {
Element facet;
facet.element = facet_insertion(f);
facet.type = type_facet;
doubleFacet(facet, doubled_nodes, doubled_facets);
}
/// double middle nodes if it's the case
if (type_facet == _segment_3) {
doubleMiddleNode(doubled_nodes, doubled_facets);
}
/// update nodal values
updateDoubledNodes(doubled_nodes);
/// loop over doubled facets to insert cohesive elements
Vector<UInt> & conn_cohesive = mesh.getConnectivity(type_cohesive);
const Vector<UInt> & conn_facet = mesh_facets.getConnectivity(type_facet);
UInt nb_nodes_per_facet = conn_facet.getNbComponent();
const Vector<Real> & position = mesh.getNodes();
Vector<UInt> & element_mat = getElementMaterial(type_cohesive);
for (UInt f = 0; f < doubled_facets.getSize(); ++f) {
UInt nb_cohesive_elements = conn_cohesive.getSize();
conn_cohesive.resize(nb_cohesive_elements + 1);
UInt first_facet = doubled_facets(f, 0);
UInt second_facet = doubled_facets(f, 1);
/// copy first facet's connectivity
for (UInt n = 0; n < nb_nodes_per_facet; ++n)
conn_cohesive(nb_cohesive_elements, n) = conn_facet(first_facet, n);
/// check if first nodes of the two facets are coincident or not
UInt first_facet_node = conn_facet(first_facet, 0);
UInt second_facet_node = conn_facet(second_facet, 0);
bool second_facet_inversion = false;
for (UInt dim = 0; dim < mesh.getSpatialDimension(); ++dim) {
if (position(first_facet_node, dim) != position(second_facet_node, dim)) {
second_facet_inversion = true;
break;
}
}
/// if the two nodes are coincident second facet connectivity is
/// normally copied, otherwise the copy is reverted
if (!second_facet_inversion) {
for (UInt n = 0; n < nb_nodes_per_facet; ++n)
conn_cohesive(nb_cohesive_elements, n + nb_nodes_per_facet)
= conn_facet(second_facet, n);
}
else {
for (UInt n = 0; n < nb_nodes_per_facet; ++n)
conn_cohesive(nb_cohesive_elements, n + nb_nodes_per_facet)
= conn_facet(second_facet, nb_nodes_per_facet - n - 1);
}
/// assign the cohesive material to the new element
element_mat.resize(nb_cohesive_elements + 1);
element_mat(nb_cohesive_elements) = cohesive_index;
materials[cohesive_index]->addElement(type_cohesive, nb_cohesive_elements, _not_ghost);
}
/// update shape functions
getFEM("CohesiveFEM").initShapeFunctions(_not_ghost);
/// resize cohesive material vectors
MaterialCohesive * mat_cohesive
= dynamic_cast<MaterialCohesive*>(materials[cohesive_index]);
AKANTU_DEBUG_ASSERT(mat_cohesive, "No cohesive materials in the materials vector");
mat_cohesive->resizeCohesiveVectors();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::doubleMiddleNode(Vector<UInt> & doubled_nodes, const Vector<UInt> & doubled_facets) {
AKANTU_DEBUG_IN();
Vector<UInt> & conn_facet = mesh_facets.getConnectivity(type_facet);
Vector<Real> & position = mesh.getNodes();
Vector<Vector<Element> > & elem_to_facet
= mesh_facets.getElementToSubelement(type_facet);
for (UInt f = 0; f < doubled_facets.getSize(); ++f) {
UInt facet_first = doubled_facets(f, 0);
UInt facet_second = doubled_facets(f, 1);
UInt new_node = position.getSize();
/// store doubled nodes
UInt nb_doubled_nodes = doubled_nodes.getSize();
doubled_nodes.resize(nb_doubled_nodes + 1);
UInt old_node = conn_facet(facet_first, 2);
doubled_nodes(nb_doubled_nodes, 0) = old_node;
doubled_nodes(nb_doubled_nodes, 1) = new_node;
/// update position
position.resize(new_node + 1);
for (UInt dim = 0; dim < spatial_dimension; ++dim)
position(new_node, dim) = position(old_node, dim);
/// update facet connectivity
conn_facet(facet_second, 2) = new_node;
/// update element connectivity
for (UInt el = 0; el < elem_to_facet(facet_second).getSize(); ++el) {
const ElementType type_elem = elem_to_facet(facet_second)(el).type;
if (type_elem != _not_defined) {
UInt elem_global = elem_to_facet(facet_second)(el).element;
Vector<UInt> & conn_elem = mesh.getConnectivity(type_elem);
for (UInt n = 0; n < conn_elem.getNbComponent(); ++n) {
if (conn_elem(elem_global, n) == old_node)
conn_elem(elem_global, n) = new_node;
}
}
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::updateDoubledNodes(const Vector<UInt> & doubled_nodes) {
AKANTU_DEBUG_IN();
UInt nb_old_nodes = displacement->getSize();
UInt nb_new_nodes = nb_old_nodes + doubled_nodes.getSize();
displacement ->resize(nb_new_nodes);
velocity ->resize(nb_new_nodes);
acceleration ->resize(nb_new_nodes);
increment_acceleration->resize(nb_new_nodes);
force ->resize(nb_new_nodes);
residual ->resize(nb_new_nodes);
boundary ->resize(nb_new_nodes);
mass ->resize(nb_new_nodes);
/**
* @todo temporary patch, should be done in a better way that works
* always (pbc, parallel, ...)
**/
delete dof_synchronizer;
dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension);
dof_synchronizer->initLocalDOFEquationNumbers();
dof_synchronizer->initGlobalDOFEquationNumbers();
for (UInt n = 0; n < doubled_nodes.getSize(); ++n) {
UInt old_node = doubled_nodes(n, 0);
UInt new_node = doubled_nodes(n, 1);
for (UInt dim = 0; dim < displacement->getNbComponent(); ++dim) {
(*displacement)(new_node, dim) = (*displacement)(old_node, dim);
}
for (UInt dim = 0; dim < velocity->getNbComponent(); ++dim) {
(*velocity)(new_node, dim) = (*velocity)(old_node, dim);
}
for (UInt dim = 0; dim < acceleration->getNbComponent(); ++dim) {
(*acceleration)(new_node, dim) = (*acceleration)(old_node, dim);
}
for (UInt dim = 0; dim < increment_acceleration->getNbComponent(); ++dim) {
(*increment_acceleration)(new_node, dim)
= (*increment_acceleration)(old_node, dim);
}
}
assembleMassLumped();
updateCurrentPosition();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::doubleFacet(Element & facet,
Vector<UInt> & doubled_nodes,
Vector<UInt> & doubled_facets) {
AKANTU_DEBUG_IN();
const UInt f_index = facet.element;
const ElementType type_facet = facet.type;
const ElementType type_subfacet = mesh.getFacetElementType(type_facet);
const UInt nb_subfacet = mesh.getNbFacetsPerElement(type_facet);
Vector<Vector<Element> > & facet_to_subfacet
= mesh_facets.getElementToSubelement(type_subfacet);
Vector<Vector<Element> > & element_to_facet
= mesh_facets.getElementToSubelement(type_facet);
Vector<Element> & subfacet_to_facet
= mesh_facets.getSubelementToElement(type_facet);
/// adding a new facet by copying original one
/// create new connectivity
Vector<UInt> & conn_facet = mesh_facets.getConnectivity(type_facet);
UInt nb_facet = conn_facet.getSize();
conn_facet.resize(nb_facet + 1);
for (UInt n = 0; n < conn_facet.getNbComponent(); ++n)
conn_facet(nb_facet, n) = conn_facet(f_index, n);
/// store doubled facets
UInt nb_doubled_facets = doubled_facets.getSize();
doubled_facets.resize(nb_doubled_facets + 1);
doubled_facets(nb_doubled_facets, 0) = f_index;
doubled_facets(nb_doubled_facets, 1) = nb_facet;
/// update elements connected to facet
element_to_facet.push_back(element_to_facet(f_index));
/// set new and original facets as boundary facets
element_to_facet(f_index)(1) = ElementNull;
element_to_facet(nb_facet)(0) = element_to_facet(nb_facet)(1);
element_to_facet(nb_facet)(1) = ElementNull;
/// update facet_to_element vector
ElementType type = element_to_facet(nb_facet)(0).type;
UInt el = element_to_facet(nb_facet)(0).element;
Vector<Element> & facet_to_element = mesh_facets.getSubelementToElement(type);
UInt i;
for (i = 0; facet_to_element(el, i).element != f_index
&& i <= facet_to_element.getNbComponent(); ++i);
facet_to_element(el, i).element = nb_facet;
/// create new links to subfacets and update list of facets
/// connected to subfacets
subfacet_to_facet.resize(nb_facet + 1);
for (UInt sf = 0; sf < nb_subfacet; ++sf) {
subfacet_to_facet(nb_facet, sf) = subfacet_to_facet(f_index, sf);
UInt sf_index = subfacet_to_facet(f_index, sf).element;
/// find index to start looping around facets connected to current
/// subfacet
UInt start = 0;
UInt nb_connected_facets = facet_to_subfacet(sf_index).getSize();
while (facet_to_subfacet(sf_index)(start).element != f_index
&& start <= facet_to_subfacet(sf_index).getSize()) ++start;
/// add the new facet to the list next to the original one
++nb_connected_facets;
facet_to_subfacet(sf_index).resize(nb_connected_facets);
for (UInt f = nb_connected_facets - 1; f > start; --f) {
facet_to_subfacet(sf_index)(f) = facet_to_subfacet(sf_index)(f - 1);
}
/// check if the new facet should be inserted before or after the
/// original one: the second element connected to both original
/// and new facet will be _not_defined, so I check if the first
/// one is equal to one of the elements connected to the following
/// facet in the facet_to_subfacet vector
UInt f_start = facet_to_subfacet(sf_index)(start).element;
UInt f_next;
if (start + 2 == nb_connected_facets)
f_next = facet_to_subfacet(sf_index)(0).element;
else
f_next = facet_to_subfacet(sf_index)(start + 2).element;
if ((element_to_facet(f_start)(0) == element_to_facet(f_next)(0))
|| ( element_to_facet(f_start)(0) == element_to_facet(f_next)(1)))
facet_to_subfacet(sf_index)(start).element = nb_facet;
else
facet_to_subfacet(sf_index)(start + 1).element = nb_facet;
/// loop on every facet connected to the current subfacet
for (UInt f = start + 2; ; ++f) {
/// reset f in order to continue looping from the beginning
if (f == nb_connected_facets) f = 0;
/// exit loop if it reaches the end
if (f == start) break;
/// if current loop facet is on the boundary, double subfacet
UInt f_global = facet_to_subfacet(sf_index)(f).element;
if (element_to_facet(f_global)(1).type == _not_defined) {
doubleSubfacet(subfacet_to_facet(f_index, sf), start, f, doubled_nodes);
break;
}
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::doubleSubfacet(const Element & subfacet,
UInt start,
UInt end,
Vector<UInt> & doubled_nodes) {
AKANTU_DEBUG_IN();
const UInt sf_index = subfacet.element;
const ElementType type_subfacet = subfacet.type;
Vector<Vector<Element> > & facet_to_subfacet
= mesh_facets.getElementToSubelement(type_subfacet);
UInt nb_subfacet = facet_to_subfacet.getSize();
Vector<UInt> & conn_point = mesh_facets.getConnectivity(_point);
Vector<Real> & position = mesh.getNodes();
/// add the new subfacet
if (spatial_dimension == 2) {
UInt new_node = position.getSize();
/// add new node in connectivity
UInt new_subfacet = conn_point.getSize();
conn_point.resize(new_subfacet + 1);
conn_point(new_subfacet) = new_node;
/// store doubled nodes
UInt nb_doubled_nodes = doubled_nodes.getSize();
doubled_nodes.resize(nb_doubled_nodes + 1);
UInt old_node = doubled_nodes(nb_doubled_nodes, 0)
= conn_point(sf_index);
doubled_nodes(nb_doubled_nodes, 1) = new_node;
/// update position
position.resize(new_node + 1);
for (UInt dim = 0; dim < spatial_dimension; ++dim)
position(new_node, dim) = position(old_node, dim);
}
/// create a vector for the new subfacet in facet_to_subfacet
facet_to_subfacet.resize(nb_subfacet + 1);
UInt nb_connected_facets = facet_to_subfacet(sf_index).getSize();
/// loop over facets from start to end
for (UInt f = start + 1; ; ++f) {
/// reset f in order to continue looping from the beginning
if (f == nb_connected_facets) f = 0;
UInt f_global = facet_to_subfacet(sf_index)(f).element;
ElementType type_facet = facet_to_subfacet(sf_index)(f).type;
Vector<UInt> & conn_facet = mesh_facets.getConnectivity(type_facet);
UInt old_node = conn_point(sf_index);
UInt new_node = conn_point(conn_point.getSize() - 1);
/// update facet connectivity
UInt i;
for (i = 0; conn_facet(f_global, i) != old_node
&& i <= conn_facet.getNbComponent(); ++i);
conn_facet(f_global, i) = new_node;
/// update element connectivity
Vector<Vector<Element> > & elem_to_facet
= mesh_facets.getElementToSubelement(type_facet);
for (UInt el = 0; el < elem_to_facet(f_global).getSize(); ++el) {
const ElementType type_elem = elem_to_facet(f_global)(el).type;
if (type_elem != _not_defined) {
UInt elem_global = elem_to_facet(f_global)(el).element;
Vector<UInt> & conn_elem = mesh.getConnectivity(type_elem);
for (UInt n = 0; n < conn_elem.getNbComponent(); ++n) {
if (conn_elem(elem_global, n) == old_node)
conn_elem(elem_global, n) = new_node;
}
}
}
/// update subfacet_to_facet vector
Vector<Element> & subfacet_to_facet
= mesh_facets.getSubelementToElement(type_facet);
for (i = 0; subfacet_to_facet(f_global, i).element != sf_index
&& i <= subfacet_to_facet.getNbComponent(); ++i);
subfacet_to_facet(f_global, i).element = nb_subfacet;
/// add current facet to facet_to_subfacet last position
Element current_facet;
current_facet.element = f_global;
current_facet.type = type_facet;
facet_to_subfacet(nb_subfacet).push_back(current_facet);
/// exit loop if it reaches the end
if (f == end) break;
}
/// rearrange the facets connected to the original subfacet and
/// compute the new number of facets connected to it
if (end < start) {
for (UInt f = 0; f < start - end; ++f)
facet_to_subfacet(sf_index)(f) = facet_to_subfacet(sf_index)(f + end + 1);
nb_connected_facets = start - end;
}
else {
for (UInt f = 1; f < nb_connected_facets - end; ++f)
facet_to_subfacet(sf_index)(start + f) = facet_to_subfacet(sf_index)(end + f);
nb_connected_facets -= end - start;
}
/// resize list of facets of the original subfacet
facet_to_subfacet(sf_index).resize(nb_connected_facets);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModelCohesive::getReversibleEnergy() {
AKANTU_DEBUG_IN();
MaterialCohesive * mat_cohesive
= dynamic_cast<MaterialCohesive*>(materials[cohesive_index]);
return mat_cohesive->getReversibleEnergy();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModelCohesive::getDissipatedEnergy() {
AKANTU_DEBUG_IN();
MaterialCohesive * mat_cohesive
= dynamic_cast<MaterialCohesive*>(materials[cohesive_index]);
return mat_cohesive->getDissipatedEnergy();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::printself(std::ostream & stream, int indent) const {
std::string space;
for(Int i = 0; i < indent; i++, space += AKANTU_INDENT);
stream << space << "SolidMechanicsModelCohesive [" << std::endl;
SolidMechanicsModel::printself(stream, indent + 1);
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
__END_AKANTU__
diff --git a/test/test_fem/test_inverse_map.cc b/test/test_fem/test_inverse_map.cc
index c54650681..9547fd526 100644
--- a/test/test_fem/test_inverse_map.cc
+++ b/test/test_fem/test_inverse_map.cc
@@ -1,103 +1,103 @@
/**
* @file test_inverse_map_XXXX.cc
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @date Mon Jul 19 10:55:49 2010
*
* @brief test of the fem class
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "fem.hh"
#include "mesh.hh"
#include "mesh_io.hh"
#include "mesh_io_msh.hh"
#include "shape_lagrange.hh"
#include "integrator_gauss.hh"
/* -------------------------------------------------------------------------- */
#include <cstdlib>
#include <fstream>
#include <iostream>
/* -------------------------------------------------------------------------- */
using namespace akantu;
int main(int argc, char *argv[]) {
akantu::initialize(argc, argv);
debug::setDebugLevel(dblTest);
const ElementType type = TYPE;
UInt dim = ElementClass<type>::getSpatialDimension();
MeshIOMSH mesh_io;
Mesh my_mesh(dim);
std::stringstream meshfilename; meshfilename << type << ".msh";
mesh_io.read(meshfilename.str(), my_mesh);
UInt nb_elements = my_mesh.getNbElement(type);
///
FEMTemplate<IntegratorGauss,ShapeLagrange> *fem =
new FEMTemplate<IntegratorGauss,ShapeLagrange>(my_mesh, dim, "my_fem");
fem->initShapeFunctions();
UInt nb_quad_points = fem->getNbQuadraturePoints(type);
/// get the quadrature points coordinates
Vector<Real> coord_on_quad(nb_quad_points*nb_elements,
my_mesh.getSpatialDimension(),
"coord_on_quad");
fem->interpolateOnQuadraturePoints(my_mesh.getNodes(),
coord_on_quad,
my_mesh.getSpatialDimension(),
type);
/// loop over the quadrature points
Vector<Real>::iterator<types::RVector> it = coord_on_quad.begin(dim);
types::RVector natural_coords(dim);
Real * quad = ElementClass<type>::getQuadraturePoints();
for(UInt el = 0 ; el < nb_elements ; ++el){
for(UInt q = 0 ; q < nb_quad_points ; ++q){
- fem->inverseMap(*it,el,type,natural_coords);
+ fem->inverseMap(*it, el, type, natural_coords);
for (UInt i = 0; i < dim; ++i) {
const Real eps = 1e-15;
AKANTU_DEBUG_ASSERT(fabs(natural_coords[i] - quad[i]) < eps,
"real coordinates inversion test failed:"
<< natural_coords[i] << " - " << quad[i]
<< " = " << natural_coords[i] - quad[i]);
}
++it;
}
}
std::cout << "inverse completed over " << nb_elements << " elements" << std::endl;
delete fem;
finalize();
return EXIT_SUCCESS;
}
diff --git a/test/test_model/test_solid_mechanics_model/test_materials/io_helper_tools.cc b/test/test_model/test_solid_mechanics_model/test_materials/io_helper_tools.cc
index 5c4f78c57..90de1b6e2 100644
--- a/test/test_model/test_solid_mechanics_model/test_materials/io_helper_tools.cc
+++ b/test/test_model/test_solid_mechanics_model/test_materials/io_helper_tools.cc
@@ -1,294 +1,293 @@
/**
* @file io_helper_tools.cc
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Fri Sep 30 11:13:01 2011
*
* @brief
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "io_helper_tools.hh"
#include "aka_common.hh"
#include "mesh.hh"
#include "solid_mechanics_model.hh"
#include "material.hh"
#include "static_communicator.hh"
// #include "dof_synchronizer.hh"
/* -------------------------------------------------------------------------- */
using namespace akantu;
/* ------------------------------------------------------------------------ */
template <ElementType type>
static iohelper::ElemType getIOHelperType() { AKANTU_DEBUG_TO_IMPLEMENT(); return iohelper::MAX_ELEM_TYPE; }
template <> iohelper::ElemType getIOHelperType<_segment_2>() { return iohelper::LINE1; }
template <> iohelper::ElemType getIOHelperType<_segment_3>() { return iohelper::LINE2; }
template <> iohelper::ElemType getIOHelperType<_triangle_3>() { return iohelper::TRIANGLE1; }
template <> iohelper::ElemType getIOHelperType<_triangle_6>() { return iohelper::TRIANGLE2; }
template <> iohelper::ElemType getIOHelperType<_quadrangle_4>() { return iohelper::QUAD1; }
template <> iohelper::ElemType getIOHelperType<_quadrangle_8>() { return iohelper::QUAD2; }
template <> iohelper::ElemType getIOHelperType<_tetrahedron_4>() { return iohelper::TETRA1; }
template <> iohelper::ElemType getIOHelperType<_tetrahedron_10>() { return iohelper::TETRA2; }
template <> iohelper::ElemType getIOHelperType<_hexahedron_8>() { return iohelper::HEX1; }
iohelper::ElemType getIOHelperType(ElementType type) {
iohelper::ElemType ioh_type = iohelper::MAX_ELEM_TYPE;
#define GET_IOHELPER_TYPE(type) \
ioh_type = getIOHelperType<type>();
AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(GET_IOHELPER_TYPE);
#undef GET_IOHELPER_TYPE
return ioh_type;
}
/* ------------------------------------------------------------------------ */
void paraviewInit(iohelper::Dumper & dumper,
const SolidMechanicsModel & model,
const ElementType & type,
const std::string & filename) {
const Mesh & mesh = model.getFEM().getMesh();
UInt spatial_dimension = mesh.getSpatialDimension(type);
UInt nb_nodes = mesh.getNbNodes();
UInt nb_element = mesh.getNbElement(type);
std::stringstream filename_sstr; filename_sstr << filename << "_" << type;
UInt whoami = StaticCommunicator::getStaticCommunicator()->whoAmI();
UInt nproc = StaticCommunicator::getStaticCommunicator()->getNbProc();
dumper.SetMode(iohelper::TEXT);
dumper.SetParallelContext(whoami, nproc);
dumper.SetPoints(mesh.getNodes().values,
spatial_dimension, nb_nodes, filename_sstr.str().c_str());
dumper.SetConnectivity((int *)mesh.getConnectivity(type).values,
getIOHelperType(type), nb_element, iohelper::C_MODE);
dumper.AddNodeDataField(model.getDisplacement().values,
spatial_dimension, "displacements");
dumper.AddNodeDataField(model.getVelocity().values,
spatial_dimension, "velocity");
dumper.AddNodeDataField(model.getAcceleration().values,
spatial_dimension, "acceleration");
dumper.AddNodeDataField(model.getMass().values,
spatial_dimension, "mass");
dumper.AddNodeDataField(model.getResidual().values,
spatial_dimension, "force");
dumper.AddNodeDataField(model.getForce().values,
spatial_dimension, "applied_force");
dumper.AddElemDataField(model.getMaterial(0).getStrain(type).values,
spatial_dimension*spatial_dimension, "strain");
dumper.AddElemDataField(model.getMaterial(0).getStress(type).values,
spatial_dimension*spatial_dimension, "stress");
-
- if(dynamic_cast<const MaterialDamage *>(&model.getMaterial(0))) {
- const MaterialDamage & mat = dynamic_cast<const MaterialDamage &>(model.getMaterial(0));
- Real * dam = mat.getDamage(type).storage();
- dumper.AddElemDataField(dam, 1, "damage");
- }
+ try {
+ const Material & mat = model.getMaterial(0);
+ const Vector<Real> & dam = mat.getVector("damage", type);
+ dumper.AddElemDataField(dam.storage(), 1, "damage");
+ } catch(...) {};
dumper.SetEmbeddedValue("displacements", 1);
dumper.SetEmbeddedValue("applied_force", 1);
dumper.SetPrefix("paraview/");
dumper.Init();
dumper.Dump();
}
/* ------------------------------------------------------------------------ */
void paraviewDump(iohelper::Dumper & dumper) {
dumper.Dump();
}
/* ------------------------------------------------------------------------ */
/* ------------------------------------------------------------------------ */
// Vector<Real> checkpoint_displacements(0, spatial_dimension);
// Vector<Real> checkpoint_velocity(0, spatial_dimension);
// Vector<Real> checkpoint_acceleration(0, spatial_dimension);
// Vector<Real> checkpoint_force(0, spatial_dimension);
// /* ------------------------------------------------------------------------ */
// void checkpointInit(iohelper::Dumper & dumper,
// const SolidMechanicsModel & model,
// const ElementType & type,
// const std::string & filename) {
// UInt whoami = StaticCommunicator::getStaticCommunicator()->whoAmI();
// UInt nproc = StaticCommunicator::getStaticCommunicator()->getNbProc();
// Vector<Real> & displacements = model.getDisplacement();
// Vector<Real> & velocity = model.getVelocity();
// Vector<Real> & acceleration = model.getAcceleration();
// Vector<Real> & force = model.getForce();
// DOFSynchronizer & dof_synchronizer = const_cast<DOFSynchronizer &>(model.getDOFSynchronizer());
// dof_synchronizer.initScatterGatherCommunicationScheme();
// if(whoami == 0){
// const Mesh & mesh = model.getFEM().getMesh();
// UInt nb_nodes = mesh.getNbGlobalNodes();
// checkpoint_displacements.resize(nb_nodes);
// checkpoint_velocity .resize(nb_nodes);
// checkpoint_acceleration .resize(nb_nodes);
// checkpoint_force .resize(nb_nodes);
// dof_synchronizer.gather(displacements, 0, &checkpoint_displacements);
// dof_synchronizer.gather(velocity , 0, &checkpoint_velocity );
// dof_synchronizer.gather(acceleration , 0, &checkpoint_acceleration );
// dof_synchronizer.gather(force , 0, &checkpoint_force );
// UInt nb_element = mesh.getNbElement(type);
// UInt spatial_dimension = mesh.getSpatialDimension(type);
// std::stringstream filename_sstr; filename_sstr << filename << "_" << type;
// dumper.SetMode(iohelper::COMPRESSED);
// dumper.SetParallelContext(whoami, nproc);
// dumper.SetPoints(mesh.getNodes().values,
// spatial_dimension, nb_nodes, filename_sstr.str().c_str());
// dumper.SetConnectivity((int *)mesh.getConnectivity(type).values,
// getIOHelperType(type), nb_element, iohelper::C_MODE);
// dumper.AddNodeDataField(checkpoint_displacements.storage(), spatial_dimension, "displacements");
// dumper.AddNodeDataField(checkpoint_velocity .storage(), spatial_dimension, "velocity");
// dumper.AddNodeDataField(checkpoint_acceleration .storage(), spatial_dimension, "acceleration");
// dumper.AddNodeDataField(checkpoint_force .storage(), spatial_dimension, "applied_force");
// dumper.SetPrefix("restart/");
// dumper.Init();
// dumper.Dump();
// } else {
// dof_synchronizer.gather(displacements, 0);
// dof_synchronizer.gather(velocity , 0);
// dof_synchronizer.gather(acceleration , 0);
// dof_synchronizer.gather(force , 0);
// }
// }
// /* ------------------------------------------------------------------------ */
// void checkpoint(iohelper::Dumper & dumper,
// const SolidMechanicsModel & model) {
// UInt whoami = StaticCommunicator::getStaticCommunicator()->whoAmI();
// DOFSynchronizer & dof_synchronizer = const_cast<DOFSynchronizer &>(model.getDOFSynchronizer());
// Vector<Real> & displacements = model.getDisplacement();
// Vector<Real> & velocity = model.getVelocity();
// Vector<Real> & acceleration = model.getAcceleration();
// Vector<Real> & force = model.getForce();
// if(whoami == 0){
// dof_synchronizer.gather(displacements, 0, &checkpoint_displacements);
// dof_synchronizer.gather(velocity , 0, &checkpoint_velocity );
// dof_synchronizer.gather(acceleration , 0, &checkpoint_acceleration );
// dof_synchronizer.gather(force , 0, &checkpoint_force );
// dumper.Dump();
// } else {
// dof_synchronizer.gather(displacements, 0);
// dof_synchronizer.gather(velocity , 0);
// dof_synchronizer.gather(acceleration , 0);
// dof_synchronizer.gather(force , 0);
// }
// }
// /* ------------------------------------------------------------------------ */
// /* ------------------------------------------------------------------------ */
// void restart(const SolidMechanicsModel & model,
// const ElementType & type,
// const std::string & filename) {
// UInt whoami = StaticCommunicator::getStaticCommunicator()->whoAmI();
// UInt nproc = StaticCommunicator::getStaticCommunicator()->getNbProc();
// DOFSynchronizer & dof_synchronizer = const_cast<DOFSynchronizer &>(model.getDOFSynchronizer());
// dof_synchronizer.initScatterGatherCommunicationScheme();
// Vector<Real> & displacements = model.getDisplacement();
// Vector<Real> & velocity = model.getVelocity();
// Vector<Real> & acceleration = model.getAcceleration();
// Vector<Real> & force = model.getForce();
// if(whoami == 0){
// const Mesh & mesh = model.getFEM().getMesh();
// UInt nb_nodes = mesh.getNbGlobalNodes();
// UInt spatial_dimension = mesh.getSpatialDimension(type);
// std::stringstream filename_sstr; filename_sstr << filename << "_" << type;
// iohelper::ReaderRestart reader;
// reader.SetMode(iohelper::COMPRESSED);
// reader.SetParallelContext(whoami, nproc);
// reader.SetPoints(filename_sstr.str().c_str());
// reader.SetConnectivity(getIOHelperType(type));
// reader.AddNodeDataField("displacements");
// reader.AddNodeDataField("velocity");
// reader.AddNodeDataField("acceleration");
// reader.AddNodeDataField("applied_force");
// reader.SetPrefix("restart/");
// reader.Init();
// reader.Read();
// Vector<Real> restart_tmp_vect(nb_nodes, spatial_dimension);
// displacements.clear();
// velocity.clear();
// acceleration.clear();
// force.clear();
// Real * tmp = reader.GetNodeDataField("displacements");
// std::copy(tmp, tmp + nb_nodes * spatial_dimension, restart_tmp_vect.storage());
// dof_synchronizer.scatter(displacements, 0, &restart_tmp_vect);
// tmp = reader.GetNodeDataField("velocity");
// std::copy(tmp, tmp + nb_nodes * spatial_dimension, restart_tmp_vect.storage());
// dof_synchronizer.scatter(velocity , 0, &restart_tmp_vect);
// tmp = reader.GetNodeDataField("acceleration");
// std::copy(tmp, tmp + nb_nodes * spatial_dimension, restart_tmp_vect.storage());
// dof_synchronizer.scatter(acceleration , 0, &restart_tmp_vect);
// tmp = reader.GetNodeDataField("applied_force");
// std::copy(tmp, tmp + nb_nodes * spatial_dimension, restart_tmp_vect.storage());
// dof_synchronizer.scatter(force , 0, &restart_tmp_vect);
// } else {
// displacements.clear();
// velocity.clear();
// acceleration.clear();
// force.clear();
// dof_synchronizer.scatter(displacements, 0);
// dof_synchronizer.scatter(velocity , 0);
// dof_synchronizer.scatter(acceleration , 0);
// dof_synchronizer.scatter(force , 0);
// }
// }
/* ------------------------------------------------------------------------ */
diff --git a/test/test_model/test_solid_mechanics_model/test_materials/local_material_damage.cc b/test/test_model/test_solid_mechanics_model/test_materials/local_material_damage.cc
index ed980b9fb..dadd1980c 100644
--- a/test/test_model/test_solid_mechanics_model/test_materials/local_material_damage.cc
+++ b/test/test_model/test_solid_mechanics_model/test_materials/local_material_damage.cc
@@ -1,154 +1,153 @@
/**
* @file local_material_damage.cc
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Marion Chambart <marion.chambart@epfl.ch>
* @date Tue Jul 27 11:53:52 2010
*
* @brief Specialization of the material class for the damage material
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "local_material_damage.hh"
#include "solid_mechanics_model.hh"
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
LocalMaterialDamage::LocalMaterialDamage(SolidMechanicsModel & model,
const ID & id) :
Material(model, id),
damage("damage", id) {
AKANTU_DEBUG_IN();
E = 0;
nu = 1./2.;
Yd = 50;
Sd = 5000;
initInternalVector(this->damage, 1);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void LocalMaterialDamage::initMaterial() {
AKANTU_DEBUG_IN();
Material::initMaterial();
resizeInternalVector(this->damage);
lambda = nu * E / ((1 + nu) * (1 - 2*nu));
mu = E / (2 * (1 + nu));
kpa = lambda + 2./3. * mu;
- is_init = true;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void LocalMaterialDamage::computeStress(ElementType el_type, GhostType ghost_type) {
AKANTU_DEBUG_IN();
Real F[3*3];
Real sigma[3*3];
resizeInternalVector(this->damage);
Real * dam = damage(el_type, ghost_type).storage();
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN;
memset(F, 0, 3 * 3 * sizeof(Real));
for (UInt i = 0; i < spatial_dimension; ++i)
for (UInt j = 0; j < spatial_dimension; ++j)
F[3*i + j] = strain_val[spatial_dimension * i + j];
for (UInt i = 0; i < spatial_dimension; ++i) F[i*3 + i] -= 1;
computeStress(F, sigma,*dam);
++dam;
for (UInt i = 0; i < spatial_dimension; ++i)
for (UInt j = 0; j < spatial_dimension; ++j)
stress_val[spatial_dimension*i + j] = sigma[3 * i + j];
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void LocalMaterialDamage::computePotentialEnergy(ElementType el_type, GhostType ghost_type) {
AKANTU_DEBUG_IN();
if(ghost_type != _not_ghost) return;
Real * epot = potential_energy(el_type).storage();
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN;
computePotentialEnergy(strain_val, stress_val, epot);
epot++;
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
bool LocalMaterialDamage::setParam(const std::string & key, const std::string & value,
const ID & id) {
std::stringstream sstr(value);
if(key == "E") { sstr >> E; }
else if(key == "nu") { sstr >> nu; }
else if(key == "Yd") { sstr >> Yd; }
else if(key == "Sd") { sstr >> Sd; }
else { return Material::setParam(key, value, id); }
return true;
}
/* -------------------------------------------------------------------------- */
void LocalMaterialDamage::printself(std::ostream & stream, int indent) const {
std::string space;
for(Int i = 0; i < indent; i++, space += AKANTU_INDENT);
stream << space << "Material<_damage> [" << std::endl;
stream << space << " + id : " << id << std::endl;
stream << space << " + name : " << name << std::endl;
stream << space << " + density : " << rho << std::endl;
stream << space << " + Young's modulus : " << E << std::endl;
stream << space << " + Poisson's ratio : " << nu << std::endl;
stream << space << " + Yd : " << Yd << std::endl;
stream << space << " + Sd : " << Sd << std::endl;
- if(is_init) {
+ if(this->isInit()) {
stream << space << " + First Lamé coefficient : " << lambda << std::endl;
stream << space << " + Second Lamé coefficient : " << mu << std::endl;
stream << space << " + Bulk coefficient : " << kpa << std::endl;
}
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
__END_AKANTU__
diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_material_damage_non_local.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_material_damage_non_local.cc
index b34300e64..f15407c51 100644
--- a/test/test_model/test_solid_mechanics_model/test_materials/test_material_damage_non_local.cc
+++ b/test/test_model/test_solid_mechanics_model/test_materials/test_material_damage_non_local.cc
@@ -1,170 +1,170 @@
/**
* @file test_material_damage_non_local.cc
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Thu Sep 8 17:31:42 2011
*
* @brief
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include <iostream>
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "mesh.hh"
#include "mesh_io.hh"
#include "mesh_io_msh.hh"
#include "solid_mechanics_model.hh"
#include "material.hh"
#include "fem.hh"
#include "local_material_damage.hh"
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_USE_IOHELPER
# include "io_helper.hh"
#endif //AKANTU_USE_IOHELPER
using namespace akantu;
akantu::Real eps = 1e-10;
static void trac(__attribute__ ((unused)) Real * position,
double * stress,
__attribute__ ((unused)) Real * normal,
__attribute__ ((unused)) UInt surface_id){
memset(stress, 0, sizeof(Real)*4);
if (fabs(position[0] - 10) < eps){
stress[0] = 3e6;
stress[3] = 3e6;
}
}
int main(int argc, char *argv[])
{
debug::setDebugLevel(dblWarning);
akantu::initialize(argc, argv);
UInt max_steps = 40000;
Real bar_height = 4.;
const UInt spatial_dimension = 2;
Mesh mesh(spatial_dimension);
MeshIOMSH mesh_io;
mesh_io.read("barre_trou.msh", mesh);
SolidMechanicsModel model(mesh);
/// model initialization
model.initVectors();
UInt nb_nodes = model.getFEM().getMesh().getNbNodes();
model.initExplicit();
model.initModel();
model.readMaterials("material_damage_non_local.dat");
model.initMaterials();
/// set vectors to 0
model.getForce().clear();
model.getVelocity().clear();
model.getAcceleration().clear();
model.getDisplacement().clear();
Real time_step = model.getStableTimeStep();
model.setTimeStep(time_step/10.);
model.assembleMassLumped();
std::cout << model << std::endl;
/// Dirichlet boundary conditions
for (akantu::UInt i = 0; i < nb_nodes; ++i) {
if(model.getFEM().getMesh().getNodes().values[spatial_dimension*i] <= eps)
model.getBoundary().values[spatial_dimension*i] = true;
if(model.getFEM().getMesh().getNodes().values[spatial_dimension*i + 1] <= eps ||
model.getFEM().getMesh().getNodes().values[spatial_dimension*i + 1] >= bar_height - eps ) {
model.getBoundary().values[spatial_dimension*i + 1] = true;
}
}
FEM & fem_boundary = model.getFEMBoundary();
fem_boundary.initShapeFunctions();
fem_boundary.computeNormalsOnControlPoints();
model.computeForcesFromFunction(trac, akantu::_bft_stress);
- MaterialDamage & mat = dynamic_cast<MaterialDamage &>((model.getMaterial(0)));
+ MaterialDamage<spatial_dimension> & mat = dynamic_cast<MaterialDamage<spatial_dimension> &>((model.getMaterial(0)));
#ifdef AKANTU_USE_IOHELPER
model.updateResidual();
iohelper::DumperParaview dumper;
dumper.SetMode(iohelper::BASE64);
dumper.SetPoints(model.getFEM().getMesh().getNodes().values, 2, nb_nodes, "coordinates_damage_nl");
dumper.SetConnectivity((int *)model.getFEM().getMesh().getConnectivity(_triangle_6).values,
iohelper::TRIANGLE2, model.getFEM().getMesh().getNbElement(_triangle_6), iohelper::C_MODE);
dumper.AddNodeDataField(model.getDisplacement().values, 2, "displacements");
dumper.AddNodeDataField(model.getVelocity().values, 2, "velocity");
dumper.AddNodeDataField(model.getForce().values, 2, "force");
dumper.AddNodeDataField(model.getMass().values, 1, "Mass");
dumper.AddNodeDataField(model.getResidual().values, 2, "residual");
dumper.AddElemDataField(mat.getStrain(_triangle_6).values, 4, "strain");
dumper.AddElemDataField(mat.getStress(_triangle_6).values, 4, "stress");
Real * dam = mat.getDamage(_triangle_6).values;
dumper.AddElemDataField(dam, 1, "damage");
dumper.SetEmbeddedValue("displacements", 1);
dumper.SetEmbeddedValue("force", 1);
dumper.SetEmbeddedValue("residual", 1);
dumper.SetEmbeddedValue("velocity", 1);
dumper.SetPrefix("paraview/");
dumper.Init();
dumper.Dump();
#endif //AKANTU_USE_IOHELPER
// mat.savePairs("cl_pairs");
for(UInt s = 0; s < max_steps; ++s) {
model.explicitPred();
model.updateResidual();
model.updateAcceleration();
model.explicitCorr();
if(s % 100 == 0) std::cout << "Step " << s+1 << "/" << max_steps <<std::endl;
#ifdef AKANTU_USE_IOHELPER
if(s % 100 == 0) dumper.Dump();
#endif //AKANTU_USE_IOHELPER
}
akantu::finalize();
return EXIT_SUCCESS;
}
diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_material_elastic_caughey_damping.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_material_elastic_caughey_damping.cc
index 29f5089fd..41b1f36b5 100644
--- a/test/test_model/test_solid_mechanics_model/test_materials/test_material_elastic_caughey_damping.cc
+++ b/test/test_model/test_solid_mechanics_model/test_materials/test_material_elastic_caughey_damping.cc
@@ -1,199 +1,199 @@
/**
* @file test_material_elastic_caughey_damping.cc
* @author David Kammer <david.kammer@epfl.ch>
* @date Fri May 27 14:36:55 2011
*
* @brief test of the material elastic caughey - physical aspecte
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include <iostream>
/* -------------------------------------------------------------------------- */
#include <limits>
#include <fstream>
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "mesh.hh"
#include "mesh_io.hh"
#include "mesh_io_msh.hh"
#include "solid_mechanics_model.hh"
#include "material.hh"
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_USE_IOHELPER
# include "io_helper_tools.hh"
#endif //AKANTU_USE_IOHELPER
using namespace akantu;
static bool testFloat(Real a, Real b, Real adm_error);
int main(int argc, char *argv[])
{
akantu::initialize(argc, argv);
akantu::debug::setDebugLevel(akantu::dblWarning);
const ElementType element_type = TYPE;
- UInt dim = Mesh::getSpatialDimension(element_type);
+ const UInt dim = ElementClass<TYPE>::getSpatialDimension();
/// load mesh
Mesh mesh(dim);
MeshIOMSH mesh_io;
std::stringstream meshname_sstr;
meshname_sstr << "single_" << element_type << ".msh";
mesh_io.read(meshname_sstr.str().c_str(), mesh);
UInt max_steps = 1000;
Real time_factor = 0.1;
UInt nb_nodes = mesh.getNbNodes();
SolidMechanicsModel model(mesh);
/* ------------------------------------------------------------------------ */
/* Initialization */
/* ------------------------------------------------------------------------ */
model.initVectors();
model.getForce().clear();
model.getVelocity().clear();
model.getAcceleration().clear();
model.getDisplacement().clear();
model.updateResidual();
model.initExplicit();
model.initModel();
model.readMaterials("material_elastic_caughey_damping.dat");
- MaterialElasticCaughey & my_mat = dynamic_cast<MaterialElasticCaughey & >(model.getMaterial(0));
+ Material & my_mat = model.getMaterial(0);
Real a_value = 1e-6;//3.836e-06;
- my_mat.setAlpha(a_value);
+ my_mat.setParam("alpha", a_value);
model.initMaterials();
std::cout << model.getMaterial(0) << std::endl;
model.assembleMassLumped();
/* ------------------------------------------------------------------------ */
/* Boundary + initial conditions */
/* ------------------------------------------------------------------------ */
Vector<UInt> the_nodes(0,1);
Real imposed_disp = 0.1;
for (UInt i = 0; i < nb_nodes; ++i) {
// block lower nodes
if(mesh.getNodes().values[i*dim+1] < 0.5) {
for (UInt j=0; j<dim; ++j)
model.getBoundary().values[dim*i + j] = true;
}
// impose displacement
else {
model.getBoundary().values[dim*i + 0] = true;
model.getDisplacement().values[dim*i + 1] = imposed_disp;
the_nodes.push_back(i);
}
}
model.updateResidual();
#ifdef AKANTU_USE_IOHELPER
iohelper::DumperParaview dumper;
paraviewInit(dumper, model, element_type, "test_mat_el_cau_dump");
#endif //AKANTU_USE_IOHELPER
/// Setting time step
Real time_step = model.getStableTimeStep() * time_factor;
std::cout << "Time Step = " << time_step << "s" << std::endl;
model.setTimeStep(time_step);
/* ------------------------------------------------------------------------ */
/* Main loop */
/* ------------------------------------------------------------------------ */
for(UInt s = 1; s <= max_steps; ++s) {
if (s%100 == 0)
std::cout << "passing step " << s << "/" << max_steps << std::endl;
model.explicitPred();
model.updateResidual();
model.updateAcceleration();
model.explicitCorr();
#ifdef AKANTU_USE_IOHELPER
if(s % 100 == 0) dumper.Dump();
#endif //AKANTU_USE_IOHELPER
}
/*
for (UInt i=0; i<the_nodes.getSize(); ++i) {
std::cout << "disp " << model.getDisplacement().values[the_nodes(i)*dim+1] << "; vel " << model.getVelocity().values[the_nodes(i)*dim+1] << std::endl;
}
*/
/* ------------------------------------------------------------------------ */
/* Test solution */
/* ------------------------------------------------------------------------ */
Real disp_tol = 1e-07;
Real velo_tol = 1e-03;
// solution triangle_3
Vector<Real> disp_triangle_3(0,1); disp_triangle_3.push_back(-0.0344941);
Vector<Real> velo_triangle_3(0,1); velo_triangle_3.push_back(-433.9);
// solution quadrangle_4
Vector<Real> disp_quadrangle_4(0,1);
disp_quadrangle_4.push_back(0.0338388);
disp_quadrangle_4.push_back(0.0338388);
Vector<Real> velo_quadrangle_4(0,1);
velo_quadrangle_4.push_back(-307.221);
velo_quadrangle_4.push_back(-307.221);
// pointer to solution
Vector<Real> * disp = NULL;
Vector<Real> * velo = NULL;
if (element_type == _triangle_3) {
disp = &disp_triangle_3;
velo = &velo_triangle_3;
}
else if (element_type == _quadrangle_4) {
disp = &disp_quadrangle_4;
velo = &velo_quadrangle_4;
}
for (UInt i=0; i<the_nodes.getSize(); ++i) {
UInt node = the_nodes.values[i];
if (!testFloat(model.getDisplacement().values[node*dim+1], disp->values[i], disp_tol)) {
std::cout << "Node " << node << " has wrong disp. Computed = " << model.getDisplacement().values[node*dim+1] << " Solution = " << disp->values[i] << std::endl;
return EXIT_FAILURE;
}
if (!testFloat(model.getVelocity().values[node*dim+1], velo->values[i], velo_tol)) {
std::cout << "Node " << node << " has wrong velo. Computed = " << model.getVelocity().values[node*dim+1] << " Solution = " << velo->values[i] << std::endl;
return EXIT_FAILURE;
}
}
finalize();
std::cout << "Patch test successful!" << std::endl;
return EXIT_SUCCESS;
}
/* -------------------------------------------------------------------------- */
bool testFloat(Real a, Real b, Real adm_error) {
if (fabs(a-b) < adm_error)
return true;
else
return false;
}
diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_material_non_local.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_material_non_local.cc
index dd59d9003..d73bd9c6d 100644
--- a/test/test_model/test_solid_mechanics_model/test_materials/test_material_non_local.cc
+++ b/test/test_model/test_solid_mechanics_model/test_materials/test_material_non_local.cc
@@ -1,185 +1,185 @@
/**
* @file test_non_local_material.cc
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Mon Aug 22 14:07:08 2011
*
* @brief test of the main part of the non local materials
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
#include <iostream>
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "mesh.hh"
#include "mesh_io.hh"
#include "solid_mechanics_model.hh"
#include "material.hh"
/* -------------------------------------------------------------------------- */
using namespace akantu;
#ifdef AKANTU_USE_IOHELPER
# include "io_helper.hh"
static void paraviewInit(iohelper::Dumper & dumper, const SolidMechanicsModel & model);
//static void paraviewDump(iohelper::Dumper & dumper);
#endif
ByElementTypeReal quadrature_points_volumes("quadrature_points_volumes", "test");
const ElementType TYPE = _triangle_6;
int main(int argc, char *argv[]) {
akantu::initialize(argc, argv);
debug::setDebugLevel(akantu::dblWarning);
- UInt spatial_dimension = 2;
+ const UInt spatial_dimension = 2;
Mesh mesh(spatial_dimension);
MeshIOMSH mesh_io;
mesh_io.read("mesh.msh", mesh);
SolidMechanicsModel model(mesh);
model.initModel();
model.initVectors();
model.readMaterials("material_non_local.dat");
model.initMaterials();
// model.getFEM().getMesh().initByElementTypeVector(quadrature_points_volumes, 1, 0);
- const MaterialNonLocal<BaseWeightFunction> & mat =
- dynamic_cast<const MaterialNonLocal<BaseWeightFunction> &>(model.getMaterial(0));
+ const MaterialNonLocal<spatial_dimension, BaseWeightFunction> & mat =
+ dynamic_cast<const MaterialNonLocal<spatial_dimension, BaseWeightFunction> &>(model.getMaterial(0));
// mat.computeQuadraturePointsNeighborhoudVolumes(quadrature_points_volumes);
Real radius = mat.getRadius();
UInt nb_element = mesh.getNbElement(TYPE);
UInt nb_tot_quad = model.getFEM().getNbQuadraturePoints(TYPE) * nb_element;
Vector<Real> quads(0, spatial_dimension);
quads.resize(nb_tot_quad);
model.getFEM().interpolateOnQuadraturePoints(mesh.getNodes(),
quads, spatial_dimension,
TYPE);
Vector<Real>::iterator<types::RVector> first_quad_1 = quads.begin(spatial_dimension);
Vector<Real>::iterator<types::RVector> last_quad_1 = quads.end(spatial_dimension);
std::ofstream pout;
pout.open("bf_pairs");
UInt q1 = 0;
Real R = mat.getRadius();
for(;first_quad_1 != last_quad_1; ++first_quad_1, ++q1) {
Vector<Real>::iterator<types::RVector> first_quad_2 = quads.begin(spatial_dimension);
//Vector<Real>::iterator<types::RVector> last_quad_2 = quads.end(spatial_dimension);
UInt q2 = 0;
for(;first_quad_2 != last_quad_1; ++first_quad_2, ++q2) {
Real d = first_quad_2->distance(*first_quad_1);
if(d <= radius) {
Real alpha = (1 - d*d/(R*R));
alpha = alpha*alpha;
pout << q1 << " " << q2 << " " << alpha << std::endl;
}
}
}
pout.close();
mat.savePairs("cl_pairs");
ByElementTypeReal constant("constant_value", "test");
mesh.initByElementTypeVector(constant, 1, 0);
Mesh::type_iterator it = mesh.firstType(spatial_dimension);
Mesh::type_iterator last_type = mesh.lastType(spatial_dimension);
for(; it != last_type; ++it) {
UInt nb_quadrature_points = model.getFEM().getNbQuadraturePoints(*it);
UInt _nb_element = mesh.getNbElement(*it);
Vector<Real> & constant_vect = constant(*it);
constant_vect.resize(_nb_element * nb_quadrature_points);
std::fill_n(constant_vect.storage(), nb_quadrature_points * _nb_element, 1.);
}
ByElementTypeReal constant_avg("constant_value_avg", "test");
mesh.initByElementTypeVector(constant_avg, 1, 0);
mat.weightedAvergageOnNeighbours(constant, constant_avg, 1);
debug::setDebugLevel(akantu::dblTest);
std::cout << constant(TYPE) << std::endl;
std::cout << constant_avg(TYPE) << std::endl;
debug::setDebugLevel(akantu::dblWarning);
#ifdef AKANTU_USE_IOHELPER
iohelper::DumperParaview dumper;
paraviewInit(dumper, model);
#endif
akantu::finalize();
return EXIT_SUCCESS;
}
/* -------------------------------------------------------------------------- */
/* iohelper::Dumper vars */
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_USE_IOHELPER
/* -------------------------------------------------------------------------- */
template <ElementType type> static iohelper::ElemType paraviewType();
template <> iohelper::ElemType paraviewType<_segment_2>() { return iohelper::LINE1; }
template <> iohelper::ElemType paraviewType<_segment_3>() { return iohelper::LINE2; }
template <> iohelper::ElemType paraviewType<_triangle_3>() { return iohelper::TRIANGLE1; }
template <> iohelper::ElemType paraviewType<_triangle_6>() { return iohelper::TRIANGLE2; }
template <> iohelper::ElemType paraviewType<_quadrangle_4>() { return iohelper::QUAD1; }
template <> iohelper::ElemType paraviewType<_tetrahedron_4>() { return iohelper::TETRA1; }
template <> iohelper::ElemType paraviewType<_tetrahedron_10>() { return iohelper::TETRA2; }
template <> iohelper::ElemType paraviewType<_hexahedron_8>() { return iohelper::HEX1; }
/* -------------------------------------------------------------------------- */
void paraviewInit(iohelper::Dumper & dumper, const SolidMechanicsModel & model) {
UInt spatial_dimension = ElementClass<TYPE>::getSpatialDimension();
UInt nb_nodes = model.getFEM().getMesh().getNbNodes();
UInt nb_element = model.getFEM().getMesh().getNbElement(TYPE);
std::stringstream filename; filename << "material_non_local_" << TYPE;
dumper.SetMode(iohelper::TEXT);
dumper.SetParallelContext(StaticCommunicator::getStaticCommunicator()->whoAmI(),
StaticCommunicator::getStaticCommunicator()->getNbProc());
dumper.SetPoints(model.getFEM().getMesh().getNodes().values,
spatial_dimension, nb_nodes, filename.str().c_str());
dumper.SetConnectivity((int *)model.getFEM().getMesh().getConnectivity(TYPE).values,
paraviewType<TYPE>(), nb_element, iohelper::C_MODE);
dumper.AddElemDataField(quadrature_points_volumes(TYPE).storage(),
1, "volume");
dumper.SetPrefix("paraview/");
dumper.Init();
dumper.Dump();
}
/* -------------------------------------------------------------------------- */
// void paraviewDump(iohelper::Dumper & dumper) {
// dumper.Dump();
// }
/* -------------------------------------------------------------------------- */
#endif
diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_material_viscoelastic.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_material_viscoelastic.cc
index 1882da759..f23e2df63 100644
--- a/test/test_model/test_solid_mechanics_model/test_materials/test_material_viscoelastic.cc
+++ b/test/test_model/test_solid_mechanics_model/test_materials/test_material_viscoelastic.cc
@@ -1,185 +1,188 @@
/**
* @file test_material_viscoelastic.cc
* @author Vlad Yastrebov <vladislav.yastrebov@epfl.ch>
* @author David Kammer <david.kammer@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Feb 9 2012
*
* @brief test of the material viscoelastic
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include <limits>
#include <fstream>
#include <sstream>
#include <iostream>
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "mesh.hh"
#include "mesh_io.hh"
#include "mesh_io_msh.hh"
#include "solid_mechanics_model.hh"
#include "material.hh"
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_USE_IOHELPER
# include "io_helper_tools.hh"
#endif //AKANTU_USE_IOHELPER
using namespace akantu;
int main(int argc, char *argv[])
{
akantu::initialize(argc, argv);
akantu::debug::setDebugLevel(akantu::dblWarning);
const ElementType element_type = TYPE;
#ifdef AKANTU_USE_IOHELPER
iohelper::ElemType ioh_type = getIOHelperType(element_type);
#endif //AKANTU_USE_IOHELPER
- UInt dim = Mesh::getSpatialDimension(element_type);
+ const UInt dim = ElementClass<TYPE>::getSpatialDimension();
/// load mesh
Mesh mesh(dim);
MeshIOMSH mesh_io;
std::stringstream meshname_sstr;
meshname_sstr << "test_material_viscoelastic_" << element_type << ".msh";
mesh_io.read(meshname_sstr.str().c_str(), mesh);
UInt max_steps = 5000;
Real time_factor = 0.8;
UInt nb_nodes = mesh.getNbNodes();
SolidMechanicsModel model(mesh);
/* ------------------------------------------------------------------------ */
/* Initialization */
/* ------------------------------------------------------------------------ */
model.initVectors();
model.getForce().clear();
model.getVelocity().clear();
model.getAcceleration().clear();
model.getDisplacement().clear();
model.updateResidual();
model.initExplicit();
model.initModel();
model.readMaterials("material_viscoelastic.dat");
model.initMaterials();
std::cout << model.getMaterial(0) << std::endl;
model.assembleMassLumped();
/* ------------------------------------------------------------------------ */
/* Boundary + initial conditions */
/* ------------------------------------------------------------------------ */
Real eps = 1e-16;
for (UInt i = 0; i < nb_nodes; ++i) {
if(mesh.getNodes().values[dim*i] >= 9)
model.getDisplacement().values[dim*i+1]
= (mesh.getNodes().values[dim*i] - 9) / 100.;
if(mesh.getNodes().values[dim*i] <= eps)
model.getBoundary().values[dim*i] = true;
/*
if(mesh.getNodes().values[dim*i + 1] <= eps ||
mesh.getNodes().values[dim*i + 1] >= 1 - eps ) {
model.getBoundary().values[dim*i + 1] = true;
}
*/
}
/// dump facet and surface information to paraview
#ifdef AKANTU_USE_IOHELPER
iohelper::DumperParaview dumper;
- // paraviewInit(dumper, model, element_type, "test_viscoelastic2");
- MaterialViscoElastic & mat = dynamic_cast<MaterialViscoElastic &>((model.getMaterial(0)));
+ // paraviewInit(dumper, model, element_type, "test_viscoelastic2");
dumper.SetMode(iohelper::BASE64);
+ Material & mat = model.getMaterial(0);
std::stringstream sstr; sstr << "test_material_viscoelastic_" << TYPE;
- dumper.SetPoints(model.getFEM().getMesh().getNodes().values, 2, nb_nodes, sstr.str());
+ dumper.SetPoints(model.getFEM().getMesh().getNodes().values, dim, nb_nodes, sstr.str());
dumper.SetConnectivity((int *)model.getFEM().getMesh().getConnectivity(TYPE).values,
ioh_type, model.getFEM().getMesh().getNbElement(TYPE), iohelper::C_MODE);
dumper.AddNodeDataField(model.getDisplacement().values, 2, "displacements");
dumper.AddNodeDataField(model.getVelocity().values, 2, "velocity");
dumper.AddNodeDataField(model.getForce().values, 2, "force");
dumper.AddNodeDataField(model.getMass().values, 1, "Mass");
dumper.AddNodeDataField(model.getResidual().values, 2, "residual");
dumper.AddElemDataField(mat.getStrain(TYPE).values, 4, "strain");
dumper.AddElemDataField(mat.getStress(TYPE).values, 4, "stress");
// Real * dam = mat.getHistoryIntegral(TYPE).values;
// dumper.AddElemDataField(dam, 4, "damage");
- dumper.AddElemDataField(mat.getHistoryIntegral(TYPE).values, 4, "history");
+ try{
+ const Vector<Real> & history = mat.getVector("history_integral", TYPE);
+ dumper.AddElemDataField(history.storage(), 4, "history");
+ } catch (...) {};
dumper.SetEmbeddedValue("displacements", 1);
dumper.SetEmbeddedValue("force", 1);
dumper.SetEmbeddedValue("residual", 1);
dumper.SetEmbeddedValue("velocity", 1);
dumper.SetPrefix("paraview/");
dumper.Init();
dumper.Dump();
#endif //AKANTU_USE_IOHELPER
std::stringstream filename_sstr;
filename_sstr << "test_material_viscoelastic_" << element_type << ".out";
std::ofstream energy;
energy.open(filename_sstr.str().c_str());
energy << "id epot ekin tot" << std::endl;
/// Setting time step
Real time_step = model.getStableTimeStep() * time_factor;
std::cout << "Time Step = " << time_step << "s" << std::endl;
model.setTimeStep(time_step);
/* ------------------------------------------------------------------------ */
/* Main loop */
/* ------------------------------------------------------------------------ */
for(UInt s = 1; s <= max_steps; ++s) {
model.explicitPred();
model.updateResidual();
model.updateAcceleration();
model.explicitCorr();
Real epot = model.getPotentialEnergy();
Real ekin = model.getKineticEnergy();
if(s % 10 == 0) {
std::cerr << "passing step " << s << "/" << max_steps << std::endl;
energy << s << " " << epot << " " << ekin << " " << epot + ekin
<< std::endl;
}
#ifdef AKANTU_USE_IOHELPER
if(s % 10 == 0) paraviewDump(dumper);
#endif //AKANTU_USE_IOHELPER
}
finalize();
return EXIT_SUCCESS;
}
diff --git a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_2_exemple_1_1.cc b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_2_exemple_1_1.cc
index de2ee0cfa..a842be6b4 100644
--- a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_2_exemple_1_1.cc
+++ b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_bernoulli_beam_2_exemple_1_1.cc
@@ -1,257 +1,257 @@
/**
* @file test_structural_mechanics_model_bernoulli_beam_2_exemple_1_1.cc
* @author Fabian Barras <fabian.barras@epfl.ch>
* @date Tue May 31 19:10:26 2011
*
* @brief Computation of the analytical exemple 1.1 in the TGC vol 6
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include <limits>
#include <fstream>
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "mesh.hh"
#include "mesh_io.hh"
#include "mesh_io_msh.hh"
#include "structural_mechanics_model.hh"
#include "material.hh"
/* -------------------------------------------------------------------------- */
#define TYPE _bernoulli_beam_2
using namespace akantu;
#ifdef AKANTU_USE_IOHELPER
# include "io_helper.hh"
static void paraviewInit(iohelper::Dumper & dumper, const StructuralMechanicsModel & model);
static void paraviewDump(iohelper::Dumper & dumper);
#endif
//Linear load function
static void lin_load(double * position, double * load,
__attribute__ ((unused)) Real * normal, __attribute__ ((unused)) UInt surface_id){
memset(load,0,sizeof(Real)*3);
if (position[0]<=10){
load[1]= -6000;
}
}
int main(int argc, char *argv[]){
initialize(argc, argv);
Mesh beams(2);
debug::setDebugLevel(dblWarning);
/* -------------------------------------------------------------------------- */
// Defining the mesh
UInt nb_nodes=401;
UInt nb_nodes_1=200;
UInt nb_nodes_2=nb_nodes-nb_nodes_1 - 1;
UInt nb_element=nb_nodes-1;
Vector<Real> & nodes = const_cast<Vector<Real> &>(beams.getNodes());
nodes.resize(nb_nodes);
- beams.addConnecticityType(_bernoulli_beam_2);
+ beams.addConnectivityType(_bernoulli_beam_2);
Vector<UInt> & connectivity = const_cast<Vector<UInt> &>(beams.getConnectivity(_bernoulli_beam_2));
connectivity.resize(nb_element);
for(UInt i=0; i<nb_nodes; ++i) {
nodes(i,1)=0;
}
for (UInt i = 0; i < nb_nodes_1; ++i) {
nodes(i,0)=10.*i/((Real)nb_nodes_1);
}
nodes(nb_nodes_1,0)=10;
for (UInt i = 0; i < nb_nodes_2; ++i) {
nodes(nb_nodes_1 + i + 1,0)=10+ 8.*(i+1)/((Real)nb_nodes_2);
}
for(UInt i=0; i<nb_element; ++i) {
connectivity(i,0)=i;
connectivity(i,1)=i+1;
}
akantu::MeshIOMSH mesh_io;
mesh_io.write("b_beam_2.msh", beams);
/* -------------------------------------------------------------------------- */
// Defining the materials
// akantu::ElementType type = akantu::_bernoulli_beam_2;
akantu::StructuralMechanicsModel * model;
model = new akantu::StructuralMechanicsModel(beams);
StructuralMaterial mat1;
mat1.E=3e10;
mat1.I=0.0025;
mat1.A=0.01;
model->addMaterial(mat1);
StructuralMaterial mat2 ;
mat2.E=3e10;
mat2.I=0.00128;
mat2.A=0.01;
model->addMaterial(mat2);
/* -------------------------------------------------------------------------- */
// Defining the forces
model->initModel();
model->initVectors();
const Real M = -3600; // Momentum at 3
Vector<Real> & forces = model->getForce();
Vector<Real> & displacement = model->getDisplacement();
Vector<bool> & boundary = model->getBoundary();
const Vector<Real> & N_M = model->getStress(_bernoulli_beam_2);
Vector<UInt> & element_material = model->getElementMaterial(_bernoulli_beam_2);
forces.clear();
displacement.clear();
for (UInt i = 0; i < nb_nodes_2; ++i) {
element_material(i+nb_nodes_1)=1;
}
forces(nb_nodes-1,2) += M;
model->computeForcesFromFunction(lin_load, akantu::_bft_traction);
/* -------------------------------------------------------------------------- */
// Defining the boundary conditions
boundary(0,0) = true;
boundary(0,1) = true;
boundary(0,2) = true;
boundary(nb_nodes_1,1) = true;
boundary(nb_nodes-1,1) = true;
/* -------------------------------------------------------------------------- */
// Solve
model->initImplicitSolver();
Real error;
model->assembleStiffnessMatrix();
model->getStiffnessMatrix().saveMatrix("Kb.mtx");
UInt count = 0;
#ifdef AKANTU_USE_IOHELPER
iohelper::DumperParaview dumper;
paraviewInit(dumper, *model);
#endif
do {
if(count != 0) std::cerr << count << " - " << error << std::endl;
model->updateResidual();
model->solve();
count++;
} while (!model->testConvergenceIncrement(1e-10, error) && count < 10);
std::cerr << count << " - " << error << std::endl;
/* -------------------------------------------------------------------------- */
// Post-Processing
model->computeStressOnQuad();
model->getStiffnessMatrix().saveMatrix("Ka.mtx");
std::cout<< " d1 = " << displacement(nb_nodes_1,2) << std::endl;
std::cout<< " d2 = " << displacement(nb_nodes-1,2) << std::endl;
std::cout<< " M1 = " << N_M(0,1) << std::endl;
std::cout<< " M2 = " << N_M(2*(nb_nodes-2),1) << std::endl;
#ifdef AKANTU_USE_IOHELPER
paraviewDump(dumper);
#endif
}
/* -------------------------------------------------------------------------- */
/* iohelper::Dumper vars */
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_USE_IOHELPER
/* -------------------------------------------------------------------------- */
template <ElementType type> static iohelper::ElemType paraviewType();
template <> static iohelper::ElemType paraviewType<_segment_2>() { return iohelper::LINE1; }
template <> static iohelper::ElemType paraviewType<_segment_3>() { return iohelper::LINE2; }
template <> static iohelper::ElemType paraviewType<_triangle_3>() { return iohelper::TRIANGLE1; }
template <> static iohelper::ElemType paraviewType<_triangle_6>() { return iohelper::TRIANGLE2; }
template <> static iohelper::ElemType paraviewType<_quadrangle_4>() { return iohelper::QUAD1; }
template <> static iohelper::ElemType paraviewType<_tetrahedron_4>() { return iohelper::TETRA1; }
template <> static iohelper::ElemType paraviewType<_tetrahedron_10>() { return iohelper::TETRA2; }
template <> static iohelper::ElemType paraviewType<_hexahedron_8>() { return iohelper::HEX1; }
template <> static iohelper::ElemType paraviewType<_bernoulli_beam_2>(){ return iohelper::BEAM2; }
/* -------------------------------------------------------------------------- */
void paraviewInit(iohelper::Dumper & dumper, const StructuralMechanicsModel & model) {
UInt spatial_dimension = ElementClass<TYPE>::getSpatialDimension();
UInt nb_nodes = model.getFEM().getMesh().getNbNodes();
UInt nb_element = model.getFEM().getMesh().getNbElement(TYPE);
std::stringstream filename; filename << "beam";
dumper.SetMode(iohelper::TEXT);
dumper.SetPoints(model.getFEM().getMesh().getNodes().values,
spatial_dimension, nb_nodes, filename.str().c_str());
dumper.SetConnectivity((int *)model.getFEM().getMesh().getConnectivity(TYPE).values,
paraviewType<TYPE>(), nb_element, iohelper::C_MODE);
dumper.AddNodeDataField(model.getDisplacement().values,
3, "displacements");
dumper.AddNodeDataField(model.getForce().values,
3, "applied_force");
dumper.AddElemDataField(model.getStress(_bernoulli_beam_2).values,
2, "stress");
dumper.SetPrefix("paraview/");
dumper.Init();
dumper.Dump();
}
/* -------------------------------------------------------------------------- */
void paraviewDump(iohelper::Dumper & dumper) {
dumper.Dump();
}
/* -------------------------------------------------------------------------- */
#endif
diff --git a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_boundary_bernoulli_beam_2.cc b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_boundary_bernoulli_beam_2.cc
index 1a1da233b..1358e7b80 100644
--- a/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_boundary_bernoulli_beam_2.cc
+++ b/test/test_model/test_structural_mechanics_model/test_structural_mechanics_model_boundary_bernoulli_beam_2.cc
@@ -1,111 +1,111 @@
/**
* @file test_structural_mechanics_model_boundary_bernoulli_beam_2.cc
* @author Fabian Barras <fabian.barras@epfl.ch>
* @date Thu May 26 12:52:54 2011
*
* @brief Test the computation of linear load
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include <limits>
#include <fstream>
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "mesh.hh"
#include "mesh_io.hh"
#include "mesh_io_msh.hh"
#include "structural_mechanics_model.hh"
/* -------------------------------------------------------------------------- */
#define TYPE _bernoulli_beam_2
using namespace akantu;
static void lin_load(double * position, double * load,
__attribute__ ((unused)) Real * normal, __attribute__ ((unused)) UInt surface_id){
memset(load,0,sizeof(Real)*3);
//load[1]= (position[0])*200000000;
load[1]= 200000000;
}
int main(int argc, char *argv[]){
initialize(argc, argv);
Mesh beams(2);
debug::setDebugLevel(dblWarning);
/* -------------------------------------------------------------------------- */
// Defining the mesh
UInt nb_nodes=2;
UInt nb_element=nb_nodes-1;
UInt nb_nodes_h=2;
UInt nb_nodes_v= nb_nodes-nb_nodes_h;
Vector<Real> & nodes = const_cast<Vector<Real> &>(beams.getNodes());
nodes.resize(nb_nodes);
- beams.addConnecticityType(_bernoulli_beam_2);
+ beams.addConnectivityType(_bernoulli_beam_2);
Vector<UInt> & connectivity = const_cast<Vector<UInt> &>(beams.getConnectivity(_bernoulli_beam_2));
connectivity.resize(nb_element);
for(UInt i=0; i<nb_nodes_h; ++i) {
nodes(i,0)=2./((Real)nb_nodes_h)*i;
nodes(i,1)=0;
}
for(UInt i=nb_nodes_h; i<nb_nodes; ++i) {
nodes(i,0)=2;
nodes(i,1)=2./((Real)nb_nodes_v)*(i-nb_nodes_h);
}
for(UInt i=0; i<nb_element; ++i) {
connectivity(i,0)=i;
connectivity(i,1)=i+1;
}
akantu::MeshIOMSH mesh_io;
mesh_io.write("b_beam_2.msh", beams);
/* -------------------------------------------------------------------------- */
// Defining the forces
// akantu::ElementType type = akantu::_bernoulli_beam_2;
akantu::StructuralMechanicsModel * model;
model = new akantu::StructuralMechanicsModel(beams);
model->initModel();
model->initVectors();
Vector<Real> & forces = model->getForce();
forces.clear();
model->computeForcesFromFunction(lin_load, akantu::_bft_traction);
}