diff --git a/cmake/AkantuBuildTreeSettings.cmake.in b/cmake/AkantuBuildTreeSettings.cmake.in
index 207b86b22..fe20feff6 100644
--- a/cmake/AkantuBuildTreeSettings.cmake.in
+++ b/cmake/AkantuBuildTreeSettings.cmake.in
@@ -1,34 +1,35 @@
#===============================================================================
# @file AkantuBuildTreeSettings.cmake.in
#
# @author Nicolas Richart <nicolas.richart@epfl.ch>
#
# @date Thu Dec 01 18:00:05 2011
#
# @brief Configuration for link with build tree
#
# @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/>.
#
#===============================================================================
set(AKANTU_INCLUDE_DIR
"@AKANTU_INCLUDE_DIRS@"
)
-set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "@CMAKE_SOURCE_DIR@/cmake")
\ No newline at end of file
+set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "@CMAKE_SOURCE_DIR@/cmake")
+set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "@CMAKE_SOURCE_DIR@/cmake/Modules")
\ No newline at end of file
diff --git a/cmake/AkantuInstall.cmake b/cmake/AkantuInstall.cmake
index ef646cd30..206de39d5 100644
--- a/cmake/AkantuInstall.cmake
+++ b/cmake/AkantuInstall.cmake
@@ -1,123 +1,123 @@
#===============================================================================
# @file AkantuInstall.cmake
#
# @author Nicolas Richart <nicolas.richart@epfl.ch>
#
# @date Wed Oct 17 15:19:18 2012
#
# @brief Create the files that allows users to link with Akantu in an other
# cmake project
#
# @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/>.
#
#===============================================================================
#===============================================================================
# Config gen for external packages
#===============================================================================
configure_file(cmake/AkantuBuildTreeSettings.cmake.in "${PROJECT_BINARY_DIR}/AkantuBuildTreeSettings.cmake" @ONLY)
file(WRITE "${PROJECT_BINARY_DIR}/AkantuConfigInclude.cmake" "
#===============================================================================
# @file AkantuConfigInclude.cmake
# @author Nicolas Richart <nicolas.richart@epfl.ch>
# @date Fri Jun 11 09:46:59 2010
#
# @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
#
#===============================================================================
")
-foreach(_option ${PACKAGE_SYSTEM_PACKAGES_NAMES_LIST_ALL})
+foreach(_option ${AKANTU_PACKAGE_SYSTEM_PACKAGES_NAMES_LIST_ALL})
list(FIND AKANTU_OPTION_LIST ${_option_name} _index)
if (_index EQUAL -1)
if(NOT "${_option}" STREQUAL "CORE")
if(NOT AKANTU_${_option})
set(AKANTU_${_option} OFF)
endif()
file(APPEND "${PROJECT_BINARY_DIR}/AkantuConfigInclude.cmake" "
set(AKANTU_HAS_${_option} ${AKANTU_${_option}})")
endif()
endif()
endforeach()
file(APPEND "${PROJECT_BINARY_DIR}/AkantuConfigInclude.cmake"
"
set(AKANTU_HAS_PARTITIONER ${AKANTU_PARTITIONER})
set(AKANTU_HAS_SOLVER ${AKANTU_SOLVER})
")
foreach(_option ${AKANTU_OPTION_LIST})
package_pkg_name(${_option} _pkg_name)
file(APPEND "${PROJECT_BINARY_DIR}/AkantuConfigInclude.cmake" "
list(APPEND AKANTU_OPTION_LIST ${_option})
set(AKANTU_USE_${_option} ${AKANTU_${_option}})")
if(${_pkg_name}_LIBRARIES)
file(APPEND "${PROJECT_BINARY_DIR}/AkantuConfigInclude.cmake" "
set(${_pkg_name}_LIBRARIES ${${_pkg_name}_LIBRARIES})")
endif()
if(${_pkg_name}_INCLUDE_DIR)
file(APPEND "${PROJECT_BINARY_DIR}/AkantuConfigInclude.cmake" "
set(${_pkg_name}_INCLUDE_DIR ${${_pkg_name}_INCLUDE_DIR})
")
endif()
endforeach()
# Create the AkantuConfig.cmake and AkantuConfigVersion files
get_filename_component(CONF_REL_INCLUDE_DIR "${CMAKE_INSTALL_PREFIX}" ABSOLUTE)
configure_file(cmake/AkantuConfig.cmake.in "${PROJECT_BINARY_DIR}/AkantuConfig.cmake" @ONLY)
configure_file(cmake/AkantuConfigVersion.cmake.in "${PROJECT_BINARY_DIR}/AkantuConfigVersion.cmake" @ONLY)
configure_file(cmake/AkantuUse.cmake "${PROJECT_BINARY_DIR}/AkantuUse.cmake" COPYONLY)
# Install the export set for use with the install-tree
install(FILES ${PROJECT_BINARY_DIR}/AkantuConfig.cmake
${PROJECT_BINARY_DIR}/AkantuConfigInclude.cmake
${PROJECT_BINARY_DIR}/AkantuConfigVersion.cmake
${PROJECT_SOURCE_DIR}/cmake/AkantuUse.cmake
DESTINATION lib/akantu
COMPONENT dev)
install(FILES
${PROJECT_SOURCE_DIR}/cmake/Modules/FindIOHelper.cmake
${PROJECT_SOURCE_DIR}/cmake/Modules/FindQVIEW.cmake
${PROJECT_SOURCE_DIR}/cmake/Modules/FindMumps.cmake
${PROJECT_SOURCE_DIR}/cmake/Modules/FindScotch.cmake
${PROJECT_SOURCE_DIR}/cmake/Modules/FindGMSH.cmake
DESTINATION lib/akantu/cmake
COMPONENT dev)
diff --git a/src/common/aka_point.hh b/src/common/aka_point.hh
index 48d0b88bb..59317ba46 100644
--- a/src/common/aka_point.hh
+++ b/src/common/aka_point.hh
@@ -1,226 +1,227 @@
/**
* @file aka_point.hh
* @author David Kammer <david.kammer@epfl.ch>
* @author Alejandro Aragon <alejandro.aragon@epfl.ch>
* @date Thu Jan 12 15:43:51 2012
*
* @brief class of a geometrical point
*
* @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_AKA_POINT_HH__
#define __AKANTU_AKA_POINT_HH__
#include "aka_common.hh"
+#include <cmath>
#include <cassert>
__BEGIN_AKANTU__
const Real inf = std::numeric_limits<Real>::infinity();
template <int d, typename T = Real>
class Point {
public:
typedef T value_type;
static int dim()
{ return d; }
Point() {
for (UInt i=0; i<d; ++i)
coord_[i] = value_type();
}
explicit Point(value_type const* coordinates) {
for (UInt i=0; i<d; ++i)
coord_[i] = coordinates[i];
}
// parameter constructor
template <typename... Args>
explicit Point(const Args&... args) {
static_assert(sizeof...(Args) <= d , "*** ERROR *** Number of arguments exceeded point dension");
std::fill_n(coord_, d, value_type());
value_type coord[] = { args... };
if (sizeof...(Args) != 0)
for (size_t i=0; i<d; ++i)
coord_[i] = i < sizeof...(Args) ? coord[i] : coord[sizeof...(Args) - 1];
}
/// less operator
bool operator<(const Point & p) const {
for (int i=0; i<d; ++i)
if (coord_[i] < p[i])
return true;
return false;
}
/// bool equal operator
bool operator==(const Point & p) const {
for (int i=0; i<d; ++i)
if (coord_[i] != p[i])
return false;
return true;
}
/// standard output stream operator
friend std::ostream& operator <<(std::ostream &os, const Point &p) {
os<<"{"<<p.coord_[0];
for (int i=1; i<d; ++i)
os<<","<<p.coord_[i];
os<<"}";
return os;
}
public:
//! Get copy of coordinate
value_type operator[] (UInt index) const {
assert(index < d);
return coord_[index];
}
//! Get write access to coordinate
value_type& operator[] (UInt index) {
assert(index < d);
return coord_[index];
}
Point& operator+=(const Point& p) {
for (int i=0; i<d; ++i)
coord_[i] += p.coord_[i];
return *this;
}
Point& operator-=(const Point& p) {
for (int i=0; i<d; ++i)
coord_[i] -= p.coord_[i];
return *this;
}
template <typename S>
Point& operator*=(S s) {
for (int i=0; i<d; ++i)
coord_[i] *= s;
return *this;
}
Point& normalize()
- { return (*this)*=(1/sqrt(sq_norm())); }
+ { return (*this)*=(1/std::sqrt(sq_norm())); }
value_type sq_norm() {
value_type r = value_type();
for (int i=0; i<d; ++i)
- r += pow(coord_[i],2);
+ r += std::pow(coord_[i],2);
return r;
}
private:
Real coord_[d];
};
template <int d, typename T>
Point<d,T> operator+(const Point<d,T>& p, const Point<d,T>& q) {
Point<d,T> r(p);
return r += q;
}
template <int d, typename T>
Point<d,T> operator-(const Point<d,T>& p, const Point<d,T>& q) {
Point<d,T> r(p);
return r -= q;
}
// used operator* for dot product
template <int d, typename T>
typename Point<d,T>::value_type operator*(const Point<d,T>& p, const Point<d,T>& q) {
typename Point<d,T>::value_type r = 0;
for (int i=0; i<d; ++i)
r += p[i] * q[i];
return r;
}
template <int d, typename T>
Point<d,T> operator*(const Point<d,T>& p, typename Point<d,T>::value_type s) {
Point<d,T> r(p);
return r *= s;
}
template <int d, typename T>
Point<d,T> operator*(typename Point<d,T>::value_type s, const Point<d,T>& p) {
Point<d,T> r(p);
return r *= s;
}
// cross product
template <typename T>
Point<3,T> cross(const Point<3,T>& o, const Point<3,T>& p) {
Point<3,T> r;
for (int i=0; i<3; ++i)
r[i] = o[(i+1)%3]*p[(i+2)%3] - o[(i+2)%3]*p[(i+1)%3];
return r;
}
template <int d>
struct Bounding_volume {
typedef Point<d> point_type;
typedef typename point_type::value_type value_type;
virtual value_type measure() const = 0; // {assert(false);}
virtual Bounding_volume* combine(const Bounding_volume&) const = 0; // {assert(false);}
virtual std::ostream& print(std::ostream& os) const = 0; //{ assert(false); }
Real last_time_;
point_type velocity_;
point_type acceleration_;
friend std::ostream& operator<<(std::ostream& os, const Bounding_volume& gv)
{ return gv.print(os); }
};
__END_AKANTU__
#endif /* __AKANTU_AKA_POINT_HH__ */
diff --git a/src/model/heat_transfer/heat_transfer_model.hh b/src/model/heat_transfer/heat_transfer_model.hh
index 3fa17ed12..05e8f18da 100644
--- a/src/model/heat_transfer/heat_transfer_model.hh
+++ b/src/model/heat_transfer/heat_transfer_model.hh
@@ -1,309 +1,306 @@
/**
* @file heat_transfer_model.hh
*
* @author Rui Wang <rui.wang@epfl.ch>
* @author Srinivasa Babu Ramisetti <srinivasa.ramisetti@epfl.ch>
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Sun May 01 19:14:43 2011
*
* @brief Model of Heat Transfer
*
* @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_HEAT_TRANSFER_MODEL_HH__
#define __AKANTU_HEAT_TRANSFER_MODEL_HH__
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "aka_memory.hh"
#include "model.hh"
#include "integrator_gauss.hh"
#include "shape_lagrange.hh"
namespace akantu {
class IntegrationScheme1stOrder;
// class Solver;
// class SparseMatrix;
}
__BEGIN_AKANTU__
class HeatTransferModel : public Model, public DataAccessor {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
typedef FEMTemplate<IntegratorGauss,ShapeLagrange> MyFEMType;
HeatTransferModel(UInt spatial_dimension,
const ID & id = "heat_transfer_model",
const MemoryID & memory_id = 0) ;
HeatTransferModel(Mesh & mesh,
UInt spatial_dimension = 0,
const ID & id = "heat_transfer_model",
const MemoryID & memory_id = 0);
virtual ~HeatTransferModel() ;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// generic function to initialize everything ready for explicit dynamics
void initFull(const std::string & material_file);
/// initialize the fem object of the boundary
void initFEMBoundary(bool create_surface = true);
/// set the parameters
bool setParam(const std::string & key, const std::string & value);
/// read one material file to instantiate all the materials
void readMaterials(const std::string & filename);
/// allocate all vectors
void initVectors();
/// register the tags associated with the parallel synchronizer
void initParallel(MeshPartition * partition, DataAccessor * data_accessor=NULL);
/// initialize the model
void initModel();
/// init PBC synchronizer
void initPBC();
/// function to print the contain of the class
virtual void printself(__attribute__ ((unused)) std::ostream & stream,
__attribute__ ((unused)) int indent = 0) const {};
/* ------------------------------------------------------------------------ */
/* Methods for explicit */
/* ------------------------------------------------------------------------ */
public:
/// compute and get the stable time step
Real getStableTimeStep();
/// compute the heat flux
void updateResidual();
/// calculate the lumped capacity vector for heat transfer problem
void assembleCapacityLumped();
/// update the temperature from the temperature rate
void explicitPred();
/// update the temperature rate from the increment
void explicitCorr();
// /// initialize the heat flux
// void initializeResidual(Vector<Real> &temp);
// /// initialize temperature
// void initializeTemperature(Vector<Real> &temp);
private:
/// solve the system in temperature rate @f$C\delta \dot T = q_{n+1} - C \dot T_{n}@f$
void solveExplicitLumped();
/// compute the heat flux on ghost types
void updateResidual(const GhostType & ghost_type);
/// calculate the lumped capacity vector for heat transfer problem (w ghosttype)
void assembleCapacityLumped(const GhostType & ghost_type);
/// compute the conductivity tensor for each quadrature point in an array
void computeConductivityOnQuadPoints(const GhostType & ghost_type);
/// compute vector k \grad T for each quadrature point
void computeKgradT(const GhostType & ghost_type);
/* ------------------------------------------------------------------------ */
/* Data Accessor inherited members */
/* ------------------------------------------------------------------------ */
public:
-
- inline UInt getNbDataToPack(const Element & element,
- SynchronizationTag tag) const;
- inline UInt getNbDataToUnpack(const Element & element,
- SynchronizationTag tag) const;
+ inline UInt getNbDataForElements(const Vector<Element> & elements,
+ SynchronizationTag tag) const;
inline void packData(CommunicationBuffer & buffer,
- const Element & element,
+ const Vector<Element> & elements,
SynchronizationTag tag) const;
inline void unpackData(CommunicationBuffer & buffer,
- const Element & element,
+ const Vector<Element> & elements,
SynchronizationTag tag);
inline UInt getNbDataToPack(SynchronizationTag tag) const;
inline UInt getNbDataToUnpack(SynchronizationTag tag) const;
inline void packData(CommunicationBuffer & buffer,
const UInt index,
SynchronizationTag tag) const;
inline void unpackData(CommunicationBuffer & buffer,
const UInt index,
SynchronizationTag tag);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
inline FEM & getFEMBoundary(std::string name = "");
/// get the dimension of the system space
AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt);
/// get the current value of the time step
AKANTU_GET_MACRO(TimeStep, time_step, Real);
/// set the value of the time step
AKANTU_SET_MACRO(TimeStep, time_step, Real);
/// get the assembled heat flux
AKANTU_GET_MACRO(Residual, *residual, Vector<Real>&);
/// get the lumped capacity
AKANTU_GET_MACRO(CapacityLumped, * capacity_lumped, Vector<Real>&);
/// get the boundary vector
AKANTU_GET_MACRO(Boundary, * boundary, Vector<bool>&);
/// get the external flux vector
AKANTU_GET_MACRO(ExternalFlux, * external_flux, Vector<Real>&);
/// get the temperature gradient
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(TemperatureGradient, temperature_gradient, Real);
/// get the conductivity on q points
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ConductivityOnQpoints, conductivity_on_qpoints, Real);
/// get the conductivity on q points
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(TemperatureOnQpoints, temperature_on_qpoints, Real);
/// internal variables
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(KGradtOnQpoints, k_gradt_on_qpoints, Real);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(IntBtKgT, int_bt_k_gT, Real);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(BtKgT, bt_k_gT, Real);
/// get the temperature
AKANTU_GET_MACRO(Temperature, *temperature, Vector<Real> &);
/// get the temperature derivative
AKANTU_GET_MACRO(TemperatureRate, *temperature_rate, Vector<Real> &);
/// get the equation number Vector<Int>
AKANTU_GET_MACRO(EquationNumber, *equation_number, const Vector<Int> &);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
IntegrationScheme1stOrder * integrator;
/// time step
Real time_step;
/// temperatures array
Vector<Real> * temperature;
/// temperatures derivatives array
Vector<Real> * temperature_rate;
/// increment array (@f$\delta \dot T@f$ or @f$\delta T@f$)
Vector<Real> * increment;
/// the spatial dimension
UInt spatial_dimension;
/// the density
Real density;
/// the speed of the changing temperature
ByElementTypeReal temperature_gradient;
/// temperature field on quadrature points
ByElementTypeReal temperature_on_qpoints;
/// conductivity tensor on quadrature points
ByElementTypeReal conductivity_on_qpoints;
/// vector k \grad T on quad points
ByElementTypeReal k_gradt_on_qpoints;
/// vector \int \grad N k \grad T
ByElementTypeReal int_bt_k_gT;
/// vector \grad N k \grad T
ByElementTypeReal bt_k_gT;
//external flux vector
Vector<Real> * external_flux;
/// residuals array
Vector<Real> * residual;
/// position of a dof in the K matrix
Vector<Int> * equation_number;
//lumped vector
Vector<Real> * capacity_lumped;
/// boundary vector
Vector<bool> * boundary;
//realtime
Real time;
///capacity
Real capacity;
//conductivity matrix
Real* conductivity;
//linear variation of the conductivity (for temperature dependent conductivity)
Real conductivity_variation;
// reference temperature for the interpretation of temperature variation
Real t_ref;
//the biggest parameter of conductivity matrix
Real conductivitymax;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#if defined (AKANTU_INCLUDE_INLINE_IMPL)
# include "heat_transfer_model_inline_impl.cc"
#endif
/// standard output stream operator
inline std::ostream & operator <<(std::ostream & stream, const HeatTransferModel & _this)
{
_this.printself(stream);
return stream;
}
__END_AKANTU__
#endif /* __AKANTU_HEAT_TRANSFER_MODEL_HH__ */
diff --git a/src/model/heat_transfer/heat_transfer_model_inline_impl.cc b/src/model/heat_transfer/heat_transfer_model_inline_impl.cc
index 609779a11..004e337b7 100644
--- a/src/model/heat_transfer/heat_transfer_model_inline_impl.cc
+++ b/src/model/heat_transfer/heat_transfer_model_inline_impl.cc
@@ -1,374 +1,299 @@
/**
* @file heat_transfer_model_inline_impl.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Srinivasa Babu Ramisetti <srinivasa.ramisetti@epfl.ch>
*
* @date Sun May 01 19:14:43 2011
*
* @brief Implementation of the inline functions of the HeatTransferModel 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/>.
*
*/
/* -------------------------------------------------------------------------- */
inline FEM & HeatTransferModel::getFEMBoundary(std::string name) {
return dynamic_cast<FEM &>(getFEMClassBoundary<MyFEMType>(name));
}
/* -------------------------------------------------------------------------- */
inline UInt HeatTransferModel::getNbDataToPack(SynchronizationTag tag) const{
AKANTU_DEBUG_IN();
UInt size = 0;
UInt nb_nodes = getFEM().getMesh().getNbNodes();
switch(tag) {
case _gst_htm_temperature:
case _gst_htm_capacity: {
size += nb_nodes * sizeof(Real);
break;
}
default: {
AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
AKANTU_DEBUG_OUT();
return size;
}
/* -------------------------------------------------------------------------- */
inline UInt HeatTransferModel::getNbDataToUnpack(SynchronizationTag tag) const{
AKANTU_DEBUG_IN();
UInt size = 0;
UInt nb_nodes = getFEM().getMesh().getNbNodes();
switch(tag) {
case _gst_htm_capacity:
case _gst_htm_temperature: {
size += nb_nodes * sizeof(Real);
break;
}
default: {
AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
AKANTU_DEBUG_OUT();
return size;
}
/* -------------------------------------------------------------------------- */
inline void HeatTransferModel::packData(CommunicationBuffer & buffer,
const UInt index,
SynchronizationTag tag) const{
AKANTU_DEBUG_IN();
switch(tag) {
case _gst_htm_capacity:
buffer << (*capacity_lumped)(index);
break;
case _gst_htm_temperature: {
buffer << (*temperature)(index);
break;
}
default: {
AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
inline void HeatTransferModel::unpackData(CommunicationBuffer & buffer,
const UInt index,
SynchronizationTag tag) {
AKANTU_DEBUG_IN();
switch(tag) {
case _gst_htm_capacity: {
buffer >> (*capacity_lumped)(index);
break;
}
case _gst_htm_temperature: {
buffer >> (*temperature)(index);
break;
}
default: {
AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
-inline UInt HeatTransferModel::getNbDataToPack(const Element & element,
- SynchronizationTag tag) const {
+inline UInt HeatTransferModel::getNbDataForElements(const Vector<Element> & elements,
+ SynchronizationTag tag) const {
AKANTU_DEBUG_IN();
- UInt size = 0;
- UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type);
-
-#ifndef AKANTU_NDEBUG
- size += spatial_dimension * sizeof(Real); /// position of the barycenter of the element (only for check)
- size += spatial_dimension * nb_nodes_per_element * sizeof(Real); /// position of the nodes of the element
-#endif
-
- switch(tag) {
- case _gst_htm_capacity: {
- size += nb_nodes_per_element * sizeof(Real); // capacity vector
- break;
- }
- case _gst_htm_temperature: {
- size += nb_nodes_per_element * sizeof(Real); // temperature
- break;
- }
- case _gst_htm_gradient_temperature: {
- size += spatial_dimension * sizeof(Real); // temperature gradient
- size += nb_nodes_per_element * sizeof(Real); // nodal temperatures
- size += spatial_dimension * nb_nodes_per_element * sizeof(Real); // shape derivatives
- break;
- }
- default: {
- AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag);
- }
- }
-
- AKANTU_DEBUG_OUT();
- return size;
-}
-
-/* -------------------------------------------------------------------------- */
-inline UInt HeatTransferModel::getNbDataToUnpack(const Element & element,
- SynchronizationTag tag) const {
- AKANTU_DEBUG_IN();
UInt size = 0;
- UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type);
+ UInt nb_nodes_per_element = 0;
+ Vector<Element>::const_iterator<Element> it = elements.begin();
+ Vector<Element>::const_iterator<Element> end = elements.end();
+ for (; it != end; ++it) {
+ const Element & el = *it;
+ nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type);
+ }
#ifndef AKANTU_NDEBUG
- size += spatial_dimension * sizeof(Real); /// position of the barycenter of the element (only for check)
+ size += elements.getSize() * spatial_dimension * sizeof(Real); /// position of the barycenter of the element (only for check)
size += spatial_dimension * nb_nodes_per_element * sizeof(Real); /// position of the nodes of the element
#endif
switch(tag) {
case _gst_htm_capacity: {
size += nb_nodes_per_element * sizeof(Real); // capacity vector
break;
}
case _gst_htm_temperature: {
size += nb_nodes_per_element * sizeof(Real); // temperature
break;
}
case _gst_htm_gradient_temperature: {
size += spatial_dimension * sizeof(Real); // temperature gradient
size += nb_nodes_per_element * sizeof(Real); // nodal temperatures
size += spatial_dimension * nb_nodes_per_element * sizeof(Real); // shape derivatives
break;
}
default: {
AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
AKANTU_DEBUG_OUT();
return size;
}
/* -------------------------------------------------------------------------- */
inline void HeatTransferModel::packData(CommunicationBuffer & buffer,
- const Element & element,
+ const Vector<Element> & elements,
SynchronizationTag tag) const {
- GhostType ghost_type = _not_ghost;
- UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type);
- UInt el_offset = element.element * nb_nodes_per_element;
- UInt * conn = getFEM().getMesh().getConnectivity(element.type, ghost_type).values;
-
#ifndef AKANTU_NDEBUG
- types::RVector barycenter(spatial_dimension);
- getFEM().getMesh().getBarycenter(element.element, element.type, barycenter.storage(), ghost_type);
- buffer << barycenter;
- Real * nodes = getFEM().getMesh().getNodes().values;
- for (UInt n = 0; n < nb_nodes_per_element; ++n) {
- UInt offset_conn = conn[el_offset + n];
- for (UInt s = 0; s < spatial_dimension; ++s) {
- buffer << nodes[spatial_dimension*offset_conn+s];
- }
- }
+ Vector<Element>::const_iterator<Element> bit = elements.begin();
+ Vector<Element>::const_iterator<Element> bend = elements.end();
+ for (; bit != bend; ++bit) {
+ const Element & element = *bit;
+ types::RVector barycenter(spatial_dimension);
+ mesh.getBarycenter(element.element, element.type, barycenter.storage(), element.ghost_type);
+ buffer << barycenter;
+ }
+ // packNodalDataHelper(mesh.getNodes(), buffer, elements);
#endif
switch (tag){
case _gst_htm_capacity: {
- for (UInt n = 0; n < nb_nodes_per_element; ++n) {
- UInt offset_conn = conn[el_offset + n];
- buffer << (*capacity_lumped)(offset_conn);
- }
+ packNodalDataHelper(*capacity_lumped, buffer, elements);
break;
}
case _gst_htm_temperature: {
- for (UInt n = 0; n < nb_nodes_per_element; ++n) {
- UInt offset_conn = conn[el_offset + n];
- buffer << (*temperature)(offset_conn);
- }
+ packNodalDataHelper(*temperature, buffer, elements);
break;
}
case _gst_htm_gradient_temperature: {
- Vector<Real>::const_iterator<types::RVector> it_gtemp =
- temperature_gradient(element.type, ghost_type).begin(spatial_dimension);
-
- buffer << it_gtemp[element.element];
-
- for (UInt n = 0; n < nb_nodes_per_element; ++n) {
- UInt offset_conn = conn[el_offset + n];
- buffer << (*temperature)(offset_conn);
- }
-
- Vector<Real>::const_iterator<types::RMatrix> it_shaped =
- getFEM().getShapesDerivatives(element.type, ghost_type).begin(nb_nodes_per_element,spatial_dimension);
- buffer << it_shaped[element.element];
+ packElementalDataHelper(temperature_gradient, buffer, elements);
+ packNodalDataHelper(*temperature, buffer, elements);
+ // Vector<Real>::const_iterator<types::RMatrix> it_shaped =
+ // getFEM().getShapesDerivatives(element.type, ghost_type).begin(nb_nodes_per_element,spatial_dimension);
+ // buffer << it_shaped[element.element];
break;
}
default: {
AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
}
/* -------------------------------------------------------------------------- */
inline void HeatTransferModel::unpackData(CommunicationBuffer & buffer,
- const Element & element,
- SynchronizationTag tag) {
- GhostType ghost_type = _ghost;
- UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type);
- UInt el_offset = element.element * nb_nodes_per_element;
- UInt * conn = getFEM().getMesh().getConnectivity(element.type, ghost_type).values;
-
+ const Vector<Element> & elements,
+ SynchronizationTag tag) {
#ifndef AKANTU_NDEBUG
- types::RVector barycenter_loc(spatial_dimension);
- getFEM().getMesh().getBarycenter(element.element, element.type, barycenter_loc.storage(), ghost_type);
-
- types::RVector barycenter(spatial_dimension);
- buffer >> barycenter;
- Real tolerance = 1e-15;
- for (UInt i = 0; i < spatial_dimension; ++i) {
- if(!(std::abs(barycenter(i) - barycenter_loc(i)) <= tolerance))
- AKANTU_EXCEPTION("Unpacking an unknown value for the element : "
- << element
- << "(barycenter[" << i << "] = " << barycenter_loc(i)
- << " and buffer[" << i << "] = " << barycenter(i) << ")");
- }
+ Vector<Element>::const_iterator<Element> bit = elements.begin();
+ Vector<Element>::const_iterator<Element> bend = elements.end();
+ for (; bit != bend; ++bit) {
+ const Element & element = *bit;
+
+ types::RVector barycenter_loc(spatial_dimension);
+ mesh.getBarycenter(element.element, element.type, barycenter_loc.storage(), element.ghost_type);
- types::RVector coords(spatial_dimension);
- Real * nodes = getFEM().getMesh().getNodes().values;
- for (UInt n = 0; n < nb_nodes_per_element; ++n) {
- buffer >> coords;
- UInt offset_conn = conn[el_offset + n];
- Real * coords_local = nodes+spatial_dimension*offset_conn;
+ types::RVector barycenter(spatial_dimension);
+ buffer >> barycenter;
+ Real tolerance = 1e-15;
for (UInt i = 0; i < spatial_dimension; ++i) {
- if(!(std::abs(coords(i) - coords_local[i]) <= tolerance))
- AKANTU_EXCEPTION("Unpacking to wrong node for the element : "
- << element
- << "(coords[" << i << "] = " << coords_local[i]
- << " and buffer[" << i << "] = " << coords(i) << ")");
+ if(!(std::abs(barycenter(i) - barycenter_loc(i)) <= tolerance))
+ AKANTU_DEBUG_ERROR("Unpacking an unknown value for the element: "
+ << element
+ << "(barycenter[" << i << "] = " << barycenter_loc(i)
+ << " and buffer[" << i << "] = " << barycenter(i) << ") - tag: " << tag);
}
}
+
+ // types::RVector coords(spatial_dimension);
+ // Real * nodes = getFEM().getMesh().getNodes().values;
+ // for (UInt n = 0; n < nb_nodes_per_element; ++n) {
+ // buffer >> coords;
+ // UInt offset_conn = conn[el_offset + n];
+ // Real * coords_local = nodes+spatial_dimension*offset_conn;
+ // for (UInt i = 0; i < spatial_dimension; ++i) {
+ // if(!(std::abs(coords(i) - coords_local[i]) <= tolerance))
+ // AKANTU_EXCEPTION("Unpacking to wrong node for the element : "
+ // << element
+ // << "(coords[" << i << "] = " << coords_local[i]
+ // << " and buffer[" << i << "] = " << coords(i) << ")");
+ // }
+ // }
#endif
switch (tag){
case _gst_htm_capacity: {
- for (UInt n = 0; n < nb_nodes_per_element; ++n) {
- UInt offset_conn = conn[el_offset + n];
- buffer >> (*capacity_lumped)(offset_conn);
- }
+ unpackNodalDataHelper(*capacity_lumped, buffer, elements);
break;
}
case _gst_htm_temperature: {
- for (UInt n = 0; n < nb_nodes_per_element; ++n) {
- UInt offset_conn = conn[el_offset + n];
- Real tbuffer;
- buffer >> tbuffer;
-#ifndef AKANTU_NDEBUG
- // if (!getFEM().getMesh().isPureGhostNode(offset_conn)){
- // if (std::abs(tbuffer - (*temperature)(offset_conn)) > 1e-15){
- // AKANTU_EXCEPTION(std::scientific << std::setprecision(20)
- // << "local node is impacted with a different value computed from a distant proc"
- // << " => divergence of trajectory detected " << std::endl
- // << tbuffer << " != " << (*temperature)(offset_conn)
- // << " diff is " << tbuffer - (*temperature)(offset_conn));
- // }
- // }
-#endif
- (*temperature)(offset_conn) = tbuffer;
- }
+ unpackNodalDataHelper(*temperature, buffer, elements);
break;
}
case _gst_htm_gradient_temperature: {
- Vector<Real>::iterator<types::RVector> it_gtemp =
- temperature_gradient(element.type, ghost_type).begin(spatial_dimension);
- types::RVector gtemp(spatial_dimension);
- types::RVector temp_nodes(nb_nodes_per_element);
- types::RMatrix shaped(nb_nodes_per_element,spatial_dimension);
-
- buffer >> gtemp;
- buffer >> temp_nodes;
- buffer >> shaped;
-
- // Real tolerance = 1e-15;
- if (!Math::are_vector_equal(spatial_dimension,gtemp.storage(),it_gtemp[element.element].storage())){
- Real dist = Math::distance_3d(gtemp.storage(), it_gtemp[element.element].storage());
- debug::debugger.getOutputStream().precision(20);
- std::stringstream temperatures_str;
- temperatures_str.precision(20);
- temperatures_str << std::scientific << "temperatures are ";
- for (UInt n = 0; n < nb_nodes_per_element; ++n) {
- UInt offset_conn = conn[el_offset + n];
- temperatures_str << (*temperature)(offset_conn) << " ";
- }
- Vector<Real>::iterator<types::RMatrix> it_shaped =
- const_cast<Vector<Real> &>(getFEM().getShapesDerivatives(element.type, ghost_type))
- .begin(nb_nodes_per_element,spatial_dimension);
-
-
- AKANTU_EXCEPTION("packed gradient do not match for element " << element.element << std::endl
- << "buffer is " << gtemp << " local is " << it_gtemp[element.element]
- << " dist is " << dist << std::endl
- << temperatures_str.str() << std::endl
- << std::scientific << std::setprecision(20)
- << " distant temperatures " << temp_nodes
- << "temperature gradient size " << temperature_gradient(element.type, ghost_type).getSize()
- << " number of ghost elements " << getFEM().getMesh().getNbElement(element.type,_ghost)
- << std::scientific << std::setprecision(20)
- << " shaped " << shaped
- << std::scientific << std::setprecision(20)
- << " local shaped " << it_shaped[element.element]);
- }
+ unpackElementalDataHelper(temperature_gradient, buffer, elements);
+ unpackNodalDataHelper(*temperature, buffer, elements);
+
+ // // Real tolerance = 1e-15;
+ // if (!Math::are_vector_equal(spatial_dimension,gtemp.storage(),it_gtemp[element.element].storage())){
+ // Real dist = Math::distance_3d(gtemp.storage(), it_gtemp[element.element].storage());
+ // debug::debugger.getOutputStream().precision(20);
+ // std::stringstream temperatures_str;
+ // temperatures_str.precision(20);
+ // temperatures_str << std::scientific << "temperatures are ";
+ // for (UInt n = 0; n < nb_nodes_per_element; ++n) {
+ // UInt offset_conn = conn[el_offset + n];
+ // temperatures_str << (*temperature)(offset_conn) << " ";
+ // }
+ // Vector<Real>::iterator<types::RMatrix> it_shaped =
+ // const_cast<Vector<Real> &>(getFEM().getShapesDerivatives(element.type, ghost_type))
+ // .begin(nb_nodes_per_element,spatial_dimension);
+
+
+ // AKANTU_EXCEPTION("packed gradient do not match for element " << element.element << std::endl
+ // << "buffer is " << gtemp << " local is " << it_gtemp[element.element]
+ // << " dist is " << dist << std::endl
+ // << temperatures_str.str() << std::endl
+ // << std::scientific << std::setprecision(20)
+ // << " distant temperatures " << temp_nodes
+ // << "temperature gradient size " << temperature_gradient(element.type, ghost_type).getSize()
+ // << " number of ghost elements " << getFEM().getMesh().getNbElement(element.type,_ghost)
+ // << std::scientific << std::setprecision(20)
+ // << " shaped " << shaped
+ // << std::scientific << std::setprecision(20)
+ // << " local shaped " << it_shaped[element.element]);
+ // }
break;
}
default: {
AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
}
/* -------------------------------------------------------------------------- */
diff --git a/src/model/integration_scheme/generalized_trapezoidal.hh b/src/model/integration_scheme/generalized_trapezoidal.hh
index 60c18b7f1..25248798f 100644
--- a/src/model/integration_scheme/generalized_trapezoidal.hh
+++ b/src/model/integration_scheme/generalized_trapezoidal.hh
@@ -1,168 +1,168 @@
/**
* @file generalized_trapezoidal.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Mon Jul 04 15:07:57 2011
*
* @brief Generalized Trapezoidal Method. This implementation is taken from
* Méthodes numériques en mécanique des solides by Alain Curnier \note{ISBN:
* 2-88074-247-1}
*
* @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_GENERALIZED_TRAPEZOIDAL_HH__
#define __AKANTU_GENERALIZED_TRAPEZOIDAL_HH__
#include "integration_scheme_1st_order.hh"
__BEGIN_AKANTU__
/**
- * The two differentiate equation (thermal and cinematic) are :
+ * The two differentiate equation (thermal and kinematic) are :
* @f{eqnarray*}{
* C\dot{u}_{n+1} + Ku_{n+1} = q_{n+1}\\
* u_{n+1} = u_{n} + (1-\alpha) \Delta t \dot{u}_{n} + \alpha \Delta t \dot{u}_{n+1}
* @f}
*
* To solve it :
* Predictor :
* @f{eqnarray*}{
* u^0_{n+1} &=& u_{n} + \Delta t v_{n} \\
* \dot{u}^0_{n+1} &=& \dot{u}_{n}
* @f}
*
* Solve :
* @f[ (a C + b K^i_{n+1}) w = q_{n+1} - f^i_{n+1} - C \dot{u}^i_{n+1} @f]
*
* Corrector :
* @f{eqnarray*}{
* \dot{u}^{i+1}_{n+1} &=& \dot{u}^{i}_{n+1} + a w \\
* u^{i+1}_{n+1} &=& u^{i}_{n+1} + b w
* @f}
*
* a and b depends on the resolution method : temperature (u) or temperature rate (@f$\dot{u}@f$)
*
* For temperature : @f$ w = \delta u, a = 1 / (\alpha \Delta t) , b = 1 @f$ @n
* For temperature rate : @f$ w = \delta \dot{u}, a = 1, b = \alpha \Delta t @f$
*/
class GeneralizedTrapezoidal : public IntegrationScheme1stOrder {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
GeneralizedTrapezoidal(Real alpha) : alpha(alpha) {};
virtual ~GeneralizedTrapezoidal() {};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
virtual void integrationSchemePred(Real delta_t,
Vector<Real> & u,
Vector<Real> & u_dot,
Vector<bool> & boundary);
virtual void integrationSchemeCorrTemp(Real delta_t,
Vector<Real> & u,
Vector<Real> & u_dot,
Vector<bool> & boundary,
Vector<Real> & delta);
virtual void integrationSchemeCorrTempRate(Real delta_t,
Vector<Real> & u,
Vector<Real> & u_dot,
Vector<bool> & boundary,
Vector<Real> & delta);
public:
/// the coeffichent @f{b@f} in the description
template<IntegrationSchemeCorrectorType type>
Real getTemperatureCoefficient(Real delta_t);
/// the coeffichent @f{a@f} in the description
template<IntegrationSchemeCorrectorType type>
Real getTemperatureRateCoefficient(Real delta_t);
private:
template<IntegrationSchemeCorrectorType type>
void integrationSchemeCorr(Real delta_t,
Vector<Real> & u,
Vector<Real> & u_dot,
Vector<bool> & boundary,
Vector<Real> & delta);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
AKANTU_GET_MACRO(Alpha, alpha, Real);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
/// the @f$\alpha@f$ parameter
const Real alpha;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#if defined (AKANTU_INCLUDE_INLINE_IMPL)
# include "generalized_trapezoidal_inline_impl.cc"
#endif
/**
* Forward Euler (explicit) -> condition on delta_t
*/
class ForwardEuler : public GeneralizedTrapezoidal {
public:
ForwardEuler() : GeneralizedTrapezoidal(0.) {};
};
/**
* Trapezoidal rule (implicit), midpoint rule or Crank-Nicolson
*/
class TrapezoidalRule1 : public GeneralizedTrapezoidal {
public:
TrapezoidalRule1() : GeneralizedTrapezoidal(.5) {};
};
/**
* Backward Euler (implicit)
*/
class BackwardEuler : public GeneralizedTrapezoidal {
public:
BackwardEuler() : GeneralizedTrapezoidal(1.) {};
};
__END_AKANTU__
#endif /* __AKANTU_GENERALIZED_TRAPEZOIDAL_HH__ */
diff --git a/src/model/model.hh b/src/model/model.hh
index 84ec1c4dd..66e93b96f 100644
--- a/src/model/model.hh
+++ b/src/model/model.hh
@@ -1,198 +1,224 @@
/**
* @file model.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author David Simon Kammer <david.kammer@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Tue Jul 27 18:15:37 2010
*
* @brief Interface of a model
*
* @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_MODEL_HH__
#define __AKANTU_MODEL_HH__
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "aka_memory.hh"
#include "mesh.hh"
#include "fem.hh"
#include "mesh_utils.hh"
#include "synchronizer_registry.hh"
#include "distributed_synchronizer.hh"
#include "static_communicator.hh"
#include "mesh_partition.hh"
#include "dof_synchronizer.hh"
#include "pbc_synchronizer.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
class Model : public Memory {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
typedef Mesh mesh_type;
Model(Mesh& mesh, const ID & id = "model",
const MemoryID & memory_id = 0);
virtual ~Model();
typedef std::map<std::string, FEM *> FEMMap;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
virtual void initModel() = 0;
/// create the synchronizer registry object
void createSynchronizerRegistry(DataAccessor * data_accessor);
/// create a parallel synchronizer and distribute the mesh
Synchronizer & createParallelSynch(MeshPartition * partition,
DataAccessor * data_accessor);
/// change local equation number so that PBC is assembled properly
void changeLocalEquationNumberforPBC(std::map<UInt,UInt> & pbc_pair,UInt dimension);
/// function to print the containt of the class
virtual void printself(std::ostream & stream, int indent = 0) const = 0;
/// initialize the model for PBC
void setPBC(UInt x, UInt y, UInt z);
void setPBC(SurfacePairList & surface_pairs,
ElementType surface_e_type);
virtual void initPBC();
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// get id of model
AKANTU_GET_MACRO(ID, id, const ID)
/// get the number of surfaces
AKANTU_GET_MACRO(Mesh, mesh, Mesh&);
/// return the object hadling equation numbers
AKANTU_GET_MACRO(DOFSynchronizer, *dof_synchronizer, const DOFSynchronizer &)
/// synchronize the boundary in case of parallel run
virtual void synchronizeBoundaries() {};
/// return the fem object associated with a provided name
inline FEM & getFEM(std::string name = "") const;
/// return the fem boundary object associated with a provided name
virtual FEM & getFEMBoundary(std::string name = "");
/// register a fem object associated with name
template <typename FEMClass> inline void registerFEMObject(const std::string & name,
Mesh & mesh,
UInt spatial_dimension);
/// unregister a fem object associated with name
inline void unRegisterFEMObject(const std::string & name);
/// return the synchronizer registry
SynchronizerRegistry & getSynchronizerRegistry();
public:
/// return the fem object associated with a provided name
template <typename FEMClass>
inline FEMClass & getFEMClass(std::string name = "") const;
/// return the fem boundary object associated with a provided name
template <typename FEMClass>
inline FEMClass & getFEMClassBoundary(std::string name = "");
/// get the pbc pairs
std::map<UInt,UInt> & getPBCPairs(){return pbc_pair;};
protected:
/// returns if node is slave in pbc
inline bool getIsPBCSlaveNode(const UInt node);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
+ template<typename T>
+ inline void packElementalDataHelper(const ByElementTypeVector<T> & data_to_pack,
+ CommunicationBuffer & buffer,
+ const Vector<Element> & elements) const;
+ template<typename T>
+ inline void unpackElementalDataHelper(ByElementTypeVector<T> & data_to_unpack,
+ CommunicationBuffer & buffer,
+ const Vector<Element> & elements) const;
+ template<typename T, bool pack_helper>
+ inline void packUnpackElementalDataHelper(ByElementTypeVector<T> & data_to_pack,
+ CommunicationBuffer & buffer,
+ const Vector<Element> & element) const;
+
+ /* -------------------------------------------------------------------------- */
+ template<typename T>
+ inline void packNodalDataHelper(Vector<T> & data_to_pack,
+ CommunicationBuffer & buffer,
+ const Vector<Element> & element) const;
+ template<typename T>
+ inline void unpackNodalDataHelper(Vector<T> & data_to_unpack,
+ CommunicationBuffer & buffer,
+ const Vector<Element> & element) const;
+ template<typename T, bool pack_helper>
+ inline void packUnpackNodalDataHelper(Vector<T> & data,
+ CommunicationBuffer & buffer,
+ const Vector<Element> & elements) const;
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// Mesh
Mesh & mesh;
/// id
ID id;
/// the main fem object present in all models
FEMMap fems;
/// the fem object present in all models for boundaries
FEMMap fems_boundary;
/// default fem object
std::string default_fem;
/// synchronizer registry
SynchronizerRegistry * synch_registry;
/// handle the equation number things
DOFSynchronizer * dof_synchronizer;
/// pbc pairs
std::map<UInt,UInt> pbc_pair;
/// flag per node to know is pbc slave
Vector<bool> is_pbc_slave_node;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#if defined (AKANTU_INCLUDE_INLINE_IMPL)
# include "model_inline_impl.cc"
#endif
/// standard output stream operator
inline std::ostream & operator <<(std::ostream & stream, const Model & _this)
{
_this.printself(stream);
return stream;
}
__END_AKANTU__
#endif /* __AKANTU_MODEL_HH__ */
diff --git a/src/model/model_inline_impl.cc b/src/model/model_inline_impl.cc
index 5ec38b3ff..31b389be0 100644
--- a/src/model/model_inline_impl.cc
+++ b/src/model/model_inline_impl.cc
@@ -1,174 +1,279 @@
/**
* @file model_inline_impl.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Wed Aug 25 08:50:54 2010
*
* @brief inline implementation of the model 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/>.
*
*/
/* -------------------------------------------------------------------------- */
inline SynchronizerRegistry & Model::getSynchronizerRegistry(){
AKANTU_DEBUG_ASSERT(synch_registry,"synchronizer registry not initialized:"
<< " did you call createSynchronizerRegistry ?");
return *synch_registry;
}
/* -------------------------------------------------------------------------- */
template <typename FEMClass>
inline FEMClass & Model::getFEMClassBoundary(std::string name) {
AKANTU_DEBUG_IN();
if (name == "") name = default_fem;
FEMMap::const_iterator it_boun = fems_boundary.find(name);
FEMClass * tmp_fem_boundary;
if (it_boun == fems_boundary.end()){
AKANTU_DEBUG_INFO("Creating FEM boundary " << name);
FEMMap::const_iterator it = fems.find(name);
AKANTU_DEBUG_ASSERT(it != fems.end(), "The FEM " << name << " is not registered");
UInt spatial_dimension = it->second->getElementDimension();
std::stringstream sstr; sstr << id << ":fem_boundary:" << name;
tmp_fem_boundary = new FEMClass(it->second->getMesh(),
spatial_dimension-1,
sstr.str(),
memory_id);
fems_boundary[name] = tmp_fem_boundary;
} else {
tmp_fem_boundary = dynamic_cast<FEMClass *>(it_boun->second);
}
AKANTU_DEBUG_OUT();
return *tmp_fem_boundary;
}
/* -------------------------------------------------------------------------- */
template <typename FEMClass>
inline FEMClass & Model::getFEMClass(std::string name) const{
AKANTU_DEBUG_IN();
if (name == "") name = default_fem;
FEMMap::const_iterator it = fems.find(name);
AKANTU_DEBUG_ASSERT(it != fems.end(), "The FEM " << name << " is not registered");
AKANTU_DEBUG_OUT();
return dynamic_cast<FEMClass &>(*(it->second));
}
/* -------------------------------------------------------------------------- */
inline void Model::unRegisterFEMObject(const std::string & name){
FEMMap::iterator it = fems.find(name);
AKANTU_DEBUG_ASSERT(it != fems.end(), "FEM object with name "
<< name << " was not found");
delete((*it).second);
fems.erase(it);
if (!fems.empty())
default_fem = (*fems.begin()).first;
}
/* -------------------------------------------------------------------------- */
template <typename FEMClass>
inline void Model::registerFEMObject(const std::string & name,
Mesh & mesh,
UInt spatial_dimension){
if (fems.size() == 0) default_fem = name;
#ifndef AKANTU_NDEBUG
FEMMap::iterator it = fems.find(name);
AKANTU_DEBUG_ASSERT(it == fems.end(), "FEM object with name "
<< name << " was already created");
#endif
std::stringstream sstr; sstr << id << ":fem:" << name;
fems[name] = new FEMClass(mesh, spatial_dimension, sstr.str(), memory_id);
// MeshUtils::buildFacets(fems[name]->getMesh());
// std::stringstream sstr2; sstr2 << id << ":fem_boundary:" << name;
// fems_boundary[name] = new FEMClass(mesh, spatial_dimension-1, sstr2.str(), memory_id);
}
/* -------------------------------------------------------------------------- */
inline FEM & Model::getFEM(std::string name) const{
AKANTU_DEBUG_IN();
if (name == "") name = default_fem;
FEMMap::const_iterator it = fems.find(name);
AKANTU_DEBUG_ASSERT(it != fems.end(),"The FEM " << name << " is not registered");
AKANTU_DEBUG_OUT();
return *(it->second);
}
/* -------------------------------------------------------------------------- */
inline FEM & Model::getFEMBoundary(std::string name){
AKANTU_DEBUG_IN();
if (name == "") name = default_fem;
FEMMap::const_iterator it = fems_boundary.find(name);
AKANTU_DEBUG_ASSERT(it != fems_boundary.end(),
"The FEM boundary " << name << " is not registered");
AKANTU_DEBUG_ASSERT(it->second != NULL,
"The FEM boundary " << name << " was not created");
AKANTU_DEBUG_OUT();
return *(it->second);
}
/* -------------------------------------------------------------------------- */
/// @todo : should merge with a single function which handles local and global
inline void Model::changeLocalEquationNumberforPBC(std::map<UInt,UInt> & pbc_pair,
UInt dimension){
for (std::map<UInt,UInt>::iterator it = pbc_pair.begin();
it != pbc_pair.end();++it) {
Int node_master = (*it).second;
Int node_slave = (*it).first;
for (UInt i = 0; i < dimension; ++i) {
(*dof_synchronizer->getLocalDOFEquationNumbersPointer())
(node_slave*dimension+i) = dimension*node_master+i;
(*dof_synchronizer->getGlobalDOFEquationNumbersPointer())
(node_slave*dimension+i) = dimension*node_master+i;
}
}
}
/* -------------------------------------------------------------------------- */
inline bool Model::getIsPBCSlaveNode(const UInt node) {
// if no pbc is defined, is_pbc_slave_node is of size zero
if (is_pbc_slave_node.getSize() == 0)
return false;
else
return is_pbc_slave_node(node);
}
+
+
+/* -------------------------------------------------------------------------- */
+template<typename T>
+inline void Model::packElementalDataHelper(const ByElementTypeVector<T> & data_to_pack,
+ CommunicationBuffer & buffer,
+ const Vector<Element> & elements) const {
+ packUnpackElementalDataHelper<T, true>(const_cast<ByElementTypeVector<T> &>(data_to_pack),
+ buffer,
+ elements);
+}
+
+/* -------------------------------------------------------------------------- */
+template<typename T>
+inline void Model::unpackElementalDataHelper(ByElementTypeVector<T> & data_to_unpack,
+ CommunicationBuffer & buffer,
+ const Vector<Element> & elements) const {
+ packUnpackElementalDataHelper<T, false>(data_to_unpack, buffer, elements);
+}
+
+/* -------------------------------------------------------------------------- */
+template<typename T, bool pack_helper>
+inline void Model::packUnpackElementalDataHelper(ByElementTypeVector<T> & data_to_pack,
+ CommunicationBuffer & buffer,
+ const Vector<Element> & element) const {
+ ElementType current_element_type = _not_defined;
+ GhostType current_ghost_type = _casper;
+ UInt nb_quad_per_elem = 0;
+ UInt nb_component = 0;
+
+ Vector<T> * vect = NULL;
+
+ Vector<Element>::const_iterator<Element> it = element.begin();
+ Vector<Element>::const_iterator<Element> end = element.end();
+ for (; it != end; ++it) {
+ const Element & el = *it;
+ if(el.type != current_element_type || el.ghost_type != current_ghost_type) {
+ current_element_type = el.type;
+ current_ghost_type = el.ghost_type;
+ vect = &data_to_pack(el.type, el.ghost_type);
+ nb_quad_per_elem = this->getFEM().getNbQuadraturePoints(el.type, el.ghost_type);
+ nb_component = vect->getNbComponent();
+ }
+
+ types::Vector<T> data(vect->storage() + el.element * nb_component * nb_quad_per_elem,
+ nb_component * nb_quad_per_elem);
+ if(pack_helper)
+ buffer << data;
+ else
+ buffer >> data;
+ }
+}
+
+/* -------------------------------------------------------------------------- */
+template<typename T>
+inline void Model::packNodalDataHelper(Vector<T> & data_to_pack,
+ CommunicationBuffer & buffer,
+ const Vector<Element> & element) const {
+ packUnpackNodalDataHelper<T, true>(data_to_pack, buffer, element);
+}
+
+/* -------------------------------------------------------------------------- */
+template<typename T>
+inline void Model::unpackNodalDataHelper(Vector<T> & data_to_unpack,
+ CommunicationBuffer & buffer,
+ const Vector<Element> & element) const {
+ packUnpackNodalDataHelper<T, false>(data_to_unpack, buffer, element);
+}
+
+/* -------------------------------------------------------------------------- */
+template<typename T, bool pack_helper>
+inline void Model::packUnpackNodalDataHelper(Vector<T> & data,
+ CommunicationBuffer & buffer,
+ const Vector<Element> & elements) const {
+ UInt nb_component = data.getNbComponent();
+ UInt nb_nodes_per_element = 0;
+
+ ElementType current_element_type = _not_defined;
+ GhostType current_ghost_type = _casper;
+ UInt * conn = NULL;
+
+ Vector<Element>::const_iterator<Element> it = elements.begin();
+ Vector<Element>::const_iterator<Element> end = elements.end();
+ for (; it != end; ++it) {
+ const Element & el = *it;
+ if(el.type != current_element_type || el.ghost_type != current_ghost_type) {
+ current_element_type = el.type;
+ current_ghost_type = el.ghost_type;
+ conn = mesh.getConnectivity(el.type, el.ghost_type).storage();
+ nb_nodes_per_element = Mesh::getNbNodesPerElement(el.type);
+ }
+
+ UInt el_offset = el.element * nb_nodes_per_element;
+ for (UInt n = 0; n < nb_nodes_per_element; ++n) {
+ UInt offset_conn = conn[el_offset + n];
+ types::Vector<T> data_vect(data.storage() + offset_conn * nb_component,
+ nb_component);
+
+ if(pack_helper)
+ buffer << data_vect;
+ else
+ buffer >> data_vect;
+ }
+ }
+}
diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh
index 7fee3f99e..7ca82d645 100644
--- a/src/model/solid_mechanics/solid_mechanics_model.hh
+++ b/src/model/solid_mechanics/solid_mechanics_model.hh
@@ -1,593 +1,575 @@
/**
* @file solid_mechanics_model.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Tue Jul 27 18:15:37 2010
*
* @brief Model of Solid Mechanics
*
* @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_SOLID_MECHANICS_MODEL_HH__
#define __AKANTU_SOLID_MECHANICS_MODEL_HH__
/* -------------------------------------------------------------------------- */
#include <fstream>
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "model.hh"
#include "data_accessor.hh"
#include "integrator_gauss.hh"
#include "shape_lagrange.hh"
#include "aka_types.hh"
#include "integration_scheme_2nd_order.hh"
#include "solver.hh"
#include "dumpable.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
class Material;
class IntegrationScheme2ndOrder;
class Contact;
class SparseMatrix;
}
__BEGIN_AKANTU__
class SolidMechanicsModel : public Model, public DataAccessor, public MeshEventHandler, public Dumpable<DumperParaview> {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
typedef FEMTemplate<IntegratorGauss,ShapeLagrange> MyFEMType;
SolidMechanicsModel(Mesh & mesh,
UInt spatial_dimension = 0,
const ID & id = "solid_mechanics_model",
const MemoryID & memory_id = 0);
virtual ~SolidMechanicsModel();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// initialize completely the model
virtual void initFull(std::string material_file,
AnalysisMethod method = _explicit_dynamic);
/// initialize the fem object needed for boundary conditions
void initFEMBoundary(bool create_surface = true);
/// register the tags associated with the parallel synchronizer
void initParallel(MeshPartition * partition, DataAccessor * data_accessor=NULL);
/// allocate all vectors
void initVectors();
/// initialize all internal arrays for materials
void initMaterials();
/// initialize the model
virtual void initModel();
/// init PBC synchronizer
void initPBC();
/// function to print the containt of the class
virtual void printself(std::ostream & stream, int indent = 0) const;
/* ------------------------------------------------------------------------ */
/* PBC */
/* ------------------------------------------------------------------------ */
public:
/// change the equation number for proper assembly when using PBC
void changeEquationNumberforPBC(std::map<UInt,UInt> & pbc_pair);
/// synchronize Residual for output
void synchronizeResidual();
protected:
/// register PBC synchronizer
void registerPBCSynchronizer();
/* ------------------------------------------------------------------------ */
/* Explicit */
/* ------------------------------------------------------------------------ */
public:
/// initialize the stuff for the explicit scheme
void initExplicit();
/// initialize the array needed by updateResidual (residual, current_position)
void initializeUpdateResidualData();
/// update the current position vector
void updateCurrentPosition();
/// assemble the residual for the explicit scheme
void updateResidual(bool need_initialize = true);
/**
* \brief compute the acceleration from the residual
* this function is the explicit equivalent to solveDynamic in implicit
* In the case of lumped mass just divide the residual by the mass
* In the case of not lumped mass call solveDynamic<_acceleration_corrector>
*/
void updateAcceleration();
/// Solve the system @f[ A x = \alpha b @f] with A a lumped matrix
void solveLumped(Vector<Real> & x,
const Vector<Real> & A,
const Vector<Real> & b,
const Vector<bool> & boundary,
Real alpha);
/// explicit integration predictor
void explicitPred();
/// explicit integration corrector
void explicitCorr();
/* ------------------------------------------------------------------------ */
/* Implicit */
/* ------------------------------------------------------------------------ */
public:
/// initialize the solver and the jacobian_matrix (called by initImplicit)
void initSolver(SolverOptions & options = _solver_no_options);
/// initialize the stuff for the implicit solver
void initImplicit(bool dynamic = false,
SolverOptions & solver_options = _solver_no_options);
/// solve Ma = f to get the initial acceleration
void initialAcceleration();
/// assemble the stiffness matrix
void assembleStiffnessMatrix();
/// solve @f[ A\delta u = f_ext - f_int @f] in displacement
void solveDynamic();
/// solve Ku = f
void solveStatic();
/// test the convergence (norm of increment)
bool testConvergenceIncrement(Real tolerance);
bool testConvergenceIncrement(Real tolerance, Real & error);
/// test the convergence (norm of residual)
bool testConvergenceResidual(Real tolerance);
bool testConvergenceResidual(Real tolerance, Real & error);
/// implicit time integration predictor
void implicitPred();
/// implicit time integration corrector
void implicitCorr();
protected:
/// finish the computation of residual to solve in increment
void updateResidualInternal();
/// compute the support reaction and store it in force
void updateSupportReaction();
/// compute A and solve @f[ A\delta u = f_ext - f_int @f]
template<NewmarkBeta::IntegrationSchemeCorrectorType type>
void solveDynamic(Vector<Real> & increment);
/* ------------------------------------------------------------------------ */
/* Explicit/Implicit */
/* ------------------------------------------------------------------------ */
public:
/// compute the stresses
void computeStresses();
/* ------------------------------------------------------------------------ */
/* Boundaries (solid_mechanics_model_boundary.cc) */
/* ------------------------------------------------------------------------ */
public:
class SurfaceLoadFunctor {
public:
virtual void traction(__attribute__ ((unused)) const types::Vector<Real> & position,
__attribute__ ((unused)) types::Vector<Real> & traction,
__attribute__ ((unused)) const types::Vector<Real> & normal,
__attribute__ ((unused)) Surface surface_id) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
virtual void stress(__attribute__ ((unused)) const types::Vector<Real> & position,
__attribute__ ((unused)) types::RMatrix & stress,
__attribute__ ((unused)) const types::Vector<Real> & normal,
__attribute__ ((unused)) Surface surface_id) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
};
/// compute force vector from a function(x,y,z) that describe stresses
void computeForcesFromFunction(BoundaryFunction in_function,
BoundaryFunctionType function_type) __attribute__((deprecated));
template<class Functor>
void computeForcesFromFunction(Functor & functor, BoundaryFunctionType function_type);
/// integrate a force on the boundary by providing a stress tensor
void computeForcesByStressTensor(const Vector<Real> & stresses,
const ElementType & type,
const GhostType & ghost_type);
/// integrate a force on the boundary by providing a traction vector
void computeForcesByTractionVector(const Vector<Real> & tractions,
const ElementType & type,
const GhostType & ghost_type);
/// synchronize the ghost element boundaries values
void synchronizeBoundaries();
/* ------------------------------------------------------------------------ */
/* Materials (solid_mechanics_model_material.cc) */
/* ------------------------------------------------------------------------ */
public:
/// register a material in the dynamic database
Material & registerNewMaterial(const ID & mat_type,
const std::string & opt_param = "");
template <typename M>
Material & registerNewCustomMaterial(const ID & mat_type,
const std::string & opt_param = "");
/// read the material files to instantiate all the materials
void readMaterials(const std::string & filename);
/// read a custom material with a keyword and class as template
template <typename M>
UInt readCustomMaterial(const std::string & filename,
const std::string & keyword);
/// Use a UIntData in the mesh to specify the material to use per element
void setMaterialIDsFromIntData(const std::string & data_name);
/* ------------------------------------------------------------------------ */
/* Mass (solid_mechanics_model_mass.cc) */
/* ------------------------------------------------------------------------ */
public:
/// assemble the lumped mass matrix
void assembleMassLumped();
/// assemble the mass matrix
void assembleMass();
protected:
/// assemble the lumped mass matrix for local and ghost elements
void assembleMassLumped(GhostType ghost_type);
void assembleMass(GhostType ghost_type);
/// fill a vector of rho
void computeRho(Vector<Real> & rho,
ElementType type,
GhostType ghost_type);
/* ------------------------------------------------------------------------ */
/* Data Accessor inherited members */
/* ------------------------------------------------------------------------ */
public:
inline virtual UInt getNbDataForElements(const Vector<Element> & elements,
SynchronizationTag tag) const;
inline virtual void packElementData(CommunicationBuffer & buffer,
const Vector<Element> & elements,
SynchronizationTag tag) const;
inline virtual void unpackElementData(CommunicationBuffer & buffer,
const Vector<Element> & elements,
SynchronizationTag tag);
inline virtual UInt getNbDataToPack(SynchronizationTag tag) const;
inline virtual UInt getNbDataToUnpack(SynchronizationTag tag) const;
inline virtual void packData(CommunicationBuffer & buffer,
const UInt index,
SynchronizationTag tag) const;
inline virtual void unpackData(CommunicationBuffer & buffer,
const UInt index,
SynchronizationTag tag);
protected:
inline void splitElementByMaterial(const Vector<Element> & elements,
Vector<Element> * elements_per_mat) const;
-
- template<typename T>
- inline void packElementDataHelper(Vector<T> & data_to_pack,
- CommunicationBuffer & buffer,
- const Vector<Element> & element) const;
-
- /* -------------------------------------------------------------------------- */
- template<typename T>
- inline void unpackElementDataHelper(Vector<T> & data_to_unpack,
- CommunicationBuffer & buffer,
- const Vector<Element> & element) const;
-
- /* -------------------------------------------------------------------------- */
- template<typename T, bool pack_helper>
- inline void packUnpackElementDataHelper(Vector<T> & data,
- CommunicationBuffer & buffer,
- const Vector<Element> & element) const;
-
/* ------------------------------------------------------------------------ */
/* Mesh Event Handler inherited members */
/* ------------------------------------------------------------------------ */
protected:
virtual void onNodesAdded (const Vector<UInt> & nodes_list);
virtual void onNodesRemoved(const Vector<UInt> & element_list,
const Vector<UInt> & new_numbering);
virtual void onElementsAdded (const Vector<Element> & nodes_list);
virtual void onElementsRemoved(const Vector<Element> & element_list,
const ByElementTypeUInt & new_numbering);
public:
virtual void addDumpField(const std::string & field_id);
virtual void addDumpFieldVector(const std::string & field_id);
virtual void addDumpFieldTensor(const std::string & field_id);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// return the dimension of the system space
AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt);
/// get the current value of the time step
AKANTU_GET_MACRO(TimeStep, time_step, Real);
/// set the value of the time step
AKANTU_SET_MACRO(TimeStep, time_step, Real);
/// get the value of the conversion from forces/ mass to acceleration
AKANTU_GET_MACRO(F_M2A, f_m2a, Real);
/// set the value of the conversion from forces/ mass to acceleration
AKANTU_SET_MACRO(F_M2A, f_m2a, Real);
/// get the SolidMechanicsModel::displacement vector
AKANTU_GET_MACRO(Displacement, *displacement, Vector<Real> &);
/// get the SolidMechanicsModel::current_position vector \warn only consistent
/// after a call to SolidMechanicsModel::updateCurrentPosition
AKANTU_GET_MACRO(CurrentPosition, *current_position, const Vector<Real> &);
/// get the SolidMechanicsModel::increment vector \warn only consistent if
/// SolidMechanicsModel::setIncrementFlagOn has been called before
AKANTU_GET_MACRO(Increment, *increment, Vector<Real> &);
/// get the lumped SolidMechanicsModel::mass vector
AKANTU_GET_MACRO(Mass, *mass, Vector<Real> &);
/// get the SolidMechanicsModel::velocity vector
AKANTU_GET_MACRO(Velocity, *velocity, Vector<Real> &);
/// get the SolidMechanicsModel::acceleration vector, updated by
/// SolidMechanicsModel::updateAcceleration
AKANTU_GET_MACRO(Acceleration, *acceleration, Vector<Real> &);
/// get the SolidMechanicsModel::force vector (boundary forces)
AKANTU_GET_MACRO(Force, *force, Vector<Real> &);
/// get the SolidMechanicsModel::residual vector, computed by
/// SolidMechanicsModel::updateResidual
AKANTU_GET_MACRO(Residual, *residual, Vector<Real> &);
/// get the SolidMechanicsModel::boundary vector
AKANTU_GET_MACRO(Boundary, *boundary, Vector<bool> &);
/// get the SolidMechnicsModel::incrementAcceleration vector
AKANTU_GET_MACRO(IncrementAcceleration, *increment_acceleration, Vector<Real> &);
/// get the value of the SolidMechanicsModel::increment_flag
AKANTU_GET_MACRO(IncrementFlag, increment_flag, bool);
/// get the SolidMechanicsModel::element_material vector corresponding to a
/// given akantu::ElementType
// AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementMaterial, element_material, UInt);
// AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementMaterial, element_material, UInt);
//AKANTU_GET_MACRO(ElementMaterial, element_material, const ByElementTypeVector<UInt> &);
/// vectors containing local material element index for each global element index
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementIndexByMaterial, element_index_by_material, UInt);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementIndexByMaterial, element_index_by_material, UInt);
AKANTU_GET_MACRO(ElementIndexByMaterial, element_index_by_material, const ByElementTypeVector<UInt> &);
/// get a particular material
inline Material & getMaterial(UInt mat_index);
inline const Material & getMaterial(UInt mat_index) const;
/// give the number of materials
inline UInt getNbMaterials() const { return materials.size(); };
/// give the material internal index from its id
Int getInternalIndexFromID(const ID & id) const;
/// compute the stable time step
Real getStableTimeStep();
/// compute the potential energy
Real getPotentialEnergy();
/// compute the kinetic energy
Real getKineticEnergy();
/// compute the external work (for impose displacement, the velocity should be given too)
Real getExternalWork();
/// get the energies
Real getEnergy(const std::string & energy_id);
/// compute the energy for energy
Real getEnergy(const std::string & energy_id, ElementType & type, UInt index);
/// set the Contact object
AKANTU_SET_MACRO(Contact, contact, Contact *);
/**
* @brief set the SolidMechanicsModel::increment_flag to on, the activate the
* update of the SolidMechanicsModel::increment vector
*/
void setIncrementFlagOn();
/// get the stiffness matrix
AKANTU_GET_MACRO(StiffnessMatrix, *stiffness_matrix, SparseMatrix &);
/// get the mass matrix
AKANTU_GET_MACRO(MassMatrix, *mass_matrix, SparseMatrix &);
inline FEM & getFEMBoundary(std::string name = "");
/// get integrator
AKANTU_GET_MACRO(Integrator, *integrator, const IntegrationScheme2ndOrder &);
/// get access to the internal solver
AKANTU_GET_MACRO(Solver, *solver, Solver &);
protected:
/// compute the stable time step
Real getStableTimeStep(const GhostType & ghost_type);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// time step
Real time_step;
/// conversion coefficient form force/mass to acceleration
Real f_m2a;
/// displacements array
Vector<Real> * displacement;
/// lumped mass array
Vector<Real> * mass;
/// velocities array
Vector<Real> * velocity;
/// accelerations array
Vector<Real> * acceleration;
/// accelerations array
Vector<Real> * increment_acceleration;
/// forces array
Vector<Real> * force;
/// residuals array
Vector<Real> * residual;
/// boundaries array
Vector<bool> * boundary;
/// array of current position used during update residual
Vector<Real> * current_position;
/// mass matrix
SparseMatrix * mass_matrix;
/// velocity damping matrix
SparseMatrix * velocity_damping_matrix;
/// stiffness matrix
SparseMatrix * stiffness_matrix;
/// jacobian matrix @f[A = cM + dD + K@f] with @f[c = \frac{1}{\beta \Delta t^2}, d = \frac{\gamma}{\beta \Delta t} @f]
SparseMatrix * jacobian_matrix;
// /// materials of all elements
// ByElementTypeUInt element_material;
/// vectors containing local material element index for each global element index
ByElementTypeUInt element_index_by_material;
/// list of used materials
std::vector<Material *> materials;
/// integration scheme of second order used
IntegrationScheme2ndOrder * integrator;
/// increment of displacement
Vector<Real> * increment;
/// flag defining if the increment must be computed or not
bool increment_flag;
/// solver for implicit
Solver * solver;
/// object to resolve the contact
Contact * contact;
/// the spatial dimension
UInt spatial_dimension;
AnalysisMethod method;
};
__END_AKANTU__
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#include "parser.hh"
#include "material.hh"
__BEGIN_AKANTU__
#ifdef AKANTU_USE_IOHELPER
# include "dumper_iohelper_tmpl_homogenizing_field.hh"
# include "dumper_iohelper_tmpl_material_internal_field.hh"
#endif
#include "solid_mechanics_model_tmpl.hh"
#if defined (AKANTU_INCLUDE_INLINE_IMPL)
# include "solid_mechanics_model_inline_impl.cc"
#endif
/// standard output stream operator
inline std::ostream & operator <<(std::ostream & stream, const SolidMechanicsModel & _this)
{
_this.printself(stream);
return stream;
}
__END_AKANTU__
#endif /* __AKANTU_SOLID_MECHANICS_MODEL_HH__ */
diff --git a/src/model/solid_mechanics/solid_mechanics_model_inline_impl.cc b/src/model/solid_mechanics/solid_mechanics_model_inline_impl.cc
index 0a49c6f18..4dec2fd56 100644
--- a/src/model/solid_mechanics/solid_mechanics_model_inline_impl.cc
+++ b/src/model/solid_mechanics/solid_mechanics_model_inline_impl.cc
@@ -1,452 +1,399 @@
/**
* @file solid_mechanics_model_inline_impl.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Wed Aug 04 10:58:42 2010
*
* @brief Implementation of the inline functions of the SolidMechanicsModel 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/>.
*
*/
/* -------------------------------------------------------------------------- */
inline Material & SolidMechanicsModel::getMaterial(UInt mat_index) {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(mat_index < materials.size(),
"The model " << id << " has no material no "<< mat_index);
AKANTU_DEBUG_OUT();
return *materials[mat_index];
}
/* -------------------------------------------------------------------------- */
inline const Material & SolidMechanicsModel::getMaterial(UInt mat_index) const {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(mat_index < materials.size(),
"The model " << id << " has no material no "<< mat_index);
AKANTU_DEBUG_OUT();
return *materials[mat_index];
}
/* -------------------------------------------------------------------------- */
inline FEM & SolidMechanicsModel::getFEMBoundary(std::string name) {
return dynamic_cast<FEM &>(getFEMClassBoundary<MyFEMType>(name));
}
/* -------------------------------------------------------------------------- */
inline void SolidMechanicsModel::splitElementByMaterial(const Vector<Element> & elements,
Vector<Element> * elements_per_mat) const {
ElementType current_element_type = _not_defined;
GhostType current_ghost_type = _casper;
UInt * elem_mat = NULL;
Vector<Element>::const_iterator<Element> it = elements.begin();
Vector<Element>::const_iterator<Element> end = elements.end();
for (; it != end; ++it) {
const Element & el = *it;
if(el.type != current_element_type || el.ghost_type != current_ghost_type) {
current_element_type = el.type;
current_ghost_type = el.ghost_type;
elem_mat = element_index_by_material(el.type, el.ghost_type).storage();
}
elements_per_mat[elem_mat[2*el.element+1]].push_back(el);
}
}
-/* -------------------------------------------------------------------------- */
-template<typename T>
-inline void SolidMechanicsModel::packElementDataHelper(Vector<T> & data_to_pack,
- CommunicationBuffer & buffer,
- const Vector<Element> & element) const {
- packUnpackElementDataHelper<T, true>(data_to_pack, buffer, element);
-}
-
-/* -------------------------------------------------------------------------- */
-template<typename T>
-inline void SolidMechanicsModel::unpackElementDataHelper(Vector<T> & data_to_unpack,
- CommunicationBuffer & buffer,
- const Vector<Element> & element) const {
- packUnpackElementDataHelper<T, false>(data_to_unpack, buffer, element);
-}
-
-/* -------------------------------------------------------------------------- */
-template<typename T, bool pack_helper>
-inline void SolidMechanicsModel::packUnpackElementDataHelper(Vector<T> & data,
- CommunicationBuffer & buffer,
- const Vector<Element> & elements) const {
- UInt nb_component = data.getNbComponent();
- UInt nb_nodes_per_element = 0;
-
- ElementType current_element_type = _not_defined;
- GhostType current_ghost_type = _casper;
- UInt * conn = NULL;
-
- Vector<Element>::const_iterator<Element> it = elements.begin();
- Vector<Element>::const_iterator<Element> end = elements.end();
- for (; it != end; ++it) {
- const Element & el = *it;
- if(el.type != current_element_type || el.ghost_type != current_ghost_type) {
- current_element_type = el.type;
- current_ghost_type = el.ghost_type;
- conn = mesh.getConnectivity(el.type, el.ghost_type).storage();
- nb_nodes_per_element = Mesh::getNbNodesPerElement(el.type);
- }
-
- UInt el_offset = el.element * nb_nodes_per_element;
- for (UInt n = 0; n < nb_nodes_per_element; ++n) {
- UInt offset_conn = conn[el_offset + n];
- types::Vector<T> data_vect(data.storage() + offset_conn * nb_component,
- nb_component);
-
- if(pack_helper)
- buffer << data_vect;
- else
- buffer >> data_vect;
- }
- }
-}
-
/* -------------------------------------------------------------------------- */
inline UInt SolidMechanicsModel::getNbDataForElements(const Vector<Element> & elements,
SynchronizationTag tag) const {
AKANTU_DEBUG_IN();
UInt size = 0;
#ifndef AKANTU_NDEBUG
size += elements.getSize() * spatial_dimension * sizeof(Real); /// position of the barycenter of the element (only for check)
#endif
UInt nb_nodes_per_element = 0;
Vector<Element>::const_iterator<Element> it = elements.begin();
Vector<Element>::const_iterator<Element> end = elements.end();
for (; it != end; ++it) {
const Element & el = *it;
nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type);
}
switch(tag) {
case _gst_smm_mass: {
size += nb_nodes_per_element * sizeof(Real) * spatial_dimension; // mass vector
break;
}
case _gst_smm_for_strain: {
size += nb_nodes_per_element * spatial_dimension * sizeof(Real); // displacement
break;
}
case _gst_smm_boundary: {
// force, displacement, boundary
size += nb_nodes_per_element * spatial_dimension * (2 * sizeof(Real) + sizeof(bool));
break;
}
default: { }
}
Vector<Element> * elements_per_mat = new Vector<Element>[materials.size()];
this->splitElementByMaterial(elements, elements_per_mat);
for (UInt i = 0; i < materials.size(); ++i) {
size += materials[i]->getNbDataForElements(elements_per_mat[i], tag);
}
delete [] elements_per_mat;
AKANTU_DEBUG_OUT();
return size;
}
/* -------------------------------------------------------------------------- */
inline void SolidMechanicsModel::packElementData(CommunicationBuffer & buffer,
const Vector<Element> & elements,
SynchronizationTag tag) const {
AKANTU_DEBUG_IN();
#ifndef AKANTU_NDEBUG
Vector<Element>::const_iterator<Element> bit = elements.begin();
Vector<Element>::const_iterator<Element> bend = elements.end();
for (; bit != bend; ++bit) {
const Element & element = *bit;
types::RVector barycenter(spatial_dimension);
mesh.getBarycenter(element.element, element.type, barycenter.storage(), element.ghost_type);
buffer << barycenter;
}
#endif
switch(tag) {
case _gst_smm_mass: {
- packElementDataHelper(*mass, buffer, elements);
+ packNodalDataHelper(*mass, buffer, elements);
break;
}
case _gst_smm_for_strain: {
- packElementDataHelper(*displacement, buffer, elements);
+ packNodalDataHelper(*displacement, buffer, elements);
break;
}
case _gst_smm_boundary: {
- packElementDataHelper(*force, buffer, elements);
- packElementDataHelper(*velocity, buffer, elements);
- packElementDataHelper(*boundary, buffer, elements);
+ packNodalDataHelper(*force, buffer, elements);
+ packNodalDataHelper(*velocity, buffer, elements);
+ packNodalDataHelper(*boundary, buffer, elements);
break;
}
default: {
}
}
Vector<Element> * elements_per_mat = new Vector<Element>[materials.size()];
splitElementByMaterial(elements, elements_per_mat);
for (UInt i = 0; i < materials.size(); ++i) {
materials[i]->packElementData(buffer, elements_per_mat[i], tag);
}
delete [] elements_per_mat;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
inline void SolidMechanicsModel::unpackElementData(CommunicationBuffer & buffer,
const Vector<Element> & elements,
SynchronizationTag tag) {
AKANTU_DEBUG_IN();
#ifndef AKANTU_NDEBUG
Vector<Element>::const_iterator<Element> bit = elements.begin();
Vector<Element>::const_iterator<Element> bend = elements.end();
for (; bit != bend; ++bit) {
const Element & element = *bit;
types::RVector barycenter_loc(spatial_dimension);
mesh.getBarycenter(element.element, element.type, barycenter_loc.storage(), element.ghost_type);
types::RVector barycenter(spatial_dimension);
buffer >> barycenter;
Real tolerance = 1e-15;
for (UInt i = 0; i < spatial_dimension; ++i) {
if(!(std::abs(barycenter(i) - barycenter_loc(i)) <= tolerance))
AKANTU_DEBUG_ERROR("Unpacking an unknown value for the element: "
<< element
<< "(barycenter[" << i << "] = " << barycenter_loc(i)
<< " and buffer[" << i << "] = " << barycenter(i) << ") - tag: " << tag);
}
}
#endif
switch(tag) {
case _gst_smm_mass: {
- unpackElementDataHelper(*mass, buffer, elements);
+ unpackNodalDataHelper(*mass, buffer, elements);
break;
}
case _gst_smm_for_strain: {
- unpackElementDataHelper(*displacement, buffer, elements);
+ unpackNodalDataHelper(*displacement, buffer, elements);
break;
}
case _gst_smm_boundary: {
- unpackElementDataHelper(*force, buffer, elements);
- unpackElementDataHelper(*velocity, buffer, elements);
- unpackElementDataHelper(*boundary, buffer, elements);
+ unpackNodalDataHelper(*force, buffer, elements);
+ unpackNodalDataHelper(*velocity, buffer, elements);
+ unpackNodalDataHelper(*boundary, buffer, elements);
break;
}
default: {
}
}
Vector<Element> * elements_per_mat = new Vector<Element>[materials.size()];
splitElementByMaterial(elements, elements_per_mat);
for (UInt i = 0; i < materials.size(); ++i) {
materials[i]->unpackElementData(buffer, elements_per_mat[i], tag);
}
delete [] elements_per_mat;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
inline UInt SolidMechanicsModel::getNbDataToPack(SynchronizationTag tag) const {
AKANTU_DEBUG_IN();
UInt size = 0;
// UInt nb_nodes = mesh.getNbNodes();
switch(tag) {
case _gst_smm_uv: {
size += sizeof(Real) * spatial_dimension * 2;
break;
}
case _gst_smm_res: {
size += sizeof(Real) * spatial_dimension;
break;
}
case _gst_smm_mass: {
size += sizeof(Real) * spatial_dimension;
break;
}
default: {
AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
AKANTU_DEBUG_OUT();
return size;
}
/* -------------------------------------------------------------------------- */
inline UInt SolidMechanicsModel::getNbDataToUnpack(SynchronizationTag tag) const {
AKANTU_DEBUG_IN();
UInt size = 0;
// UInt nb_nodes = mesh.getNbNodes();
switch(tag) {
case _gst_smm_uv: {
size += sizeof(Real) * spatial_dimension * 2;
break;
}
case _gst_smm_res: {
size += sizeof(Real) * spatial_dimension;
break;
}
case _gst_smm_mass: {
size += sizeof(Real) * spatial_dimension;
break;
}
default: {
AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
AKANTU_DEBUG_OUT();
return size;
}
/* -------------------------------------------------------------------------- */
inline void SolidMechanicsModel::packData(CommunicationBuffer & buffer,
const UInt index,
SynchronizationTag tag) const {
AKANTU_DEBUG_IN();
switch(tag) {
case _gst_smm_uv: {
Vector<Real>::iterator<types::RVector> it_disp = displacement->begin(spatial_dimension);
Vector<Real>::iterator<types::RVector> it_velo = velocity->begin(spatial_dimension);
buffer << it_disp[index];
buffer << it_velo[index];
break;
}
case _gst_smm_res: {
Vector<Real>::iterator<types::RVector> it_res = residual->begin(spatial_dimension);
buffer << it_res[index];
break;
}
case _gst_smm_mass: {
AKANTU_DEBUG_INFO("pack mass of node " << index << " which is " << (*mass)(index,0));
Vector<Real>::iterator<types::RVector> it_mass = mass->begin(spatial_dimension);
buffer << it_mass[index];
break;
}
default: {
AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
inline void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer,
const UInt index,
SynchronizationTag tag) {
AKANTU_DEBUG_IN();
switch(tag) {
case _gst_smm_uv: {
Vector<Real>::iterator<types::RVector> it_disp = displacement->begin(spatial_dimension);
Vector<Real>::iterator<types::RVector> it_velo = velocity->begin(spatial_dimension);
buffer >> it_disp[index];
buffer >> it_velo[index];
break;
}
case _gst_smm_res: {
Vector<Real>::iterator<types::RVector> it_res = residual->begin(spatial_dimension);
buffer >> it_res[index];
break;
}
case _gst_smm_mass: {
AKANTU_DEBUG_INFO("mass of node " << index << " was " << (*mass)(index,0));
Vector<Real>::iterator<types::RVector> it_mass = mass->begin(spatial_dimension);
buffer >> it_mass[index];
AKANTU_DEBUG_INFO("mass of node " << index << " is now " << (*mass)(index,0));
break;
}
default: {
AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag);
}
}
AKANTU_DEBUG_OUT();
}
__END_AKANTU__
#include "sparse_matrix.hh"
#include "solver.hh"
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
template<NewmarkBeta::IntegrationSchemeCorrectorType type>
void SolidMechanicsModel::solveDynamic(Vector<Real> & increment) {
AKANTU_DEBUG_INFO("Solving Ma + Cv + Ku = f");
NewmarkBeta * nmb_int = dynamic_cast<NewmarkBeta *>(integrator);
Real c = nmb_int->getAccelerationCoefficient<type>(time_step);
Real d = nmb_int->getVelocityCoefficient<type>(time_step);
Real e = nmb_int->getDisplacementCoefficient<type>(time_step);
// A = c M + d C + e K
jacobian_matrix->clear();
if(type != NewmarkBeta::_acceleration_corrector)
jacobian_matrix->add(*stiffness_matrix, e);
jacobian_matrix->add(*mass_matrix, c);
mass_matrix->saveMatrix("M.mtx");
if(velocity_damping_matrix)
jacobian_matrix->add(*velocity_damping_matrix, d);
jacobian_matrix->applyBoundary(*boundary);
#ifndef AKANTU_NDEBUG
if(AKANTU_DEBUG_TEST(dblDump))
jacobian_matrix->saveMatrix("J.mtx");
#endif
jacobian_matrix->saveMatrix("J.mtx");
solver->setRHS(*residual);
// solve A w = f
solver->solve(increment);
}
/* -------------------------------------------------------------------------- */
diff --git a/test/test_model/test_heat_transfer_model/test_heat_transfer_model_square2d.cc b/test/test_model/test_heat_transfer_model/test_heat_transfer_model_square2d.cc
index 31e16dbba..e2432631f 100644
--- a/test/test_model/test_heat_transfer_model/test_heat_transfer_model_square2d.cc
+++ b/test/test_model/test_heat_transfer_model/test_heat_transfer_model_square2d.cc
@@ -1,167 +1,168 @@
/**
* @file test_heat_transfer_model_square2d.cc
*
*
* @date Sun May 01 19:14:43 2011
*
* @brief test of the class HeatTransferModel on the 3d cube
*
* @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 "mesh.hh"
#include "mesh_io.hh"
#include "mesh_io_msh.hh"
#include "heat_transfer_model.hh"
#include "pbc_synchronizer.hh"
/* -------------------------------------------------------------------------- */
#include <iostream>
#include <fstream>
-#include <string.h>
+#include <string>
using namespace std;
/* -------------------------------------------------------------------------- */
-#ifdef AKANTU_USE_IOHELPER
#include "io_helper.hh"
static void paraviewInit(akantu::HeatTransferModel & model, iohelper::Dumper & dumper);
static void paraviewDump(iohelper::Dumper & dumper);
iohelper::ElemType paraview_type = iohelper::TRIANGLE1;
-#endif //AKANTU_USE_IOHELPER
-
/* -------------------------------------------------------------------------- */
akantu::UInt spatial_dimension = 2;
akantu:: ElementType type = akantu::_triangle_3;
/* -------------------------------------------------------------------------- */
+std::string base_name;
+
int main(int argc, char *argv[])
{
akantu::debug::setDebugLevel(akantu::dblWarning);
akantu::initialize(argc, argv);
//create mesh
akantu::Mesh mesh(spatial_dimension);
akantu::MeshIOMSH mesh_io;
mesh_io.read("square_tri3.msh", mesh);
akantu::HeatTransferModel model(mesh);
//initialize everything
model.initFull("material.dat");
//assemble the lumped capacity
model.assembleCapacityLumped();
//get stable time step
akantu::Real time_step = model.getStableTimeStep()*0.8;
cout<<"time step is:" << time_step << endl;
model.setTimeStep(time_step);
//boundary conditions
const akantu::Vector<akantu::Real> & nodes = model.getFEM().getMesh().getNodes();
akantu::Vector<bool> & boundary = model.getBoundary();
akantu::Vector<akantu::Real> & temperature = model.getTemperature();
double length;
length = 1.;
akantu::UInt nb_nodes = model.getFEM().getMesh().getNbNodes();
for (akantu::UInt i = 0; i < nb_nodes; ++i) {
temperature(i) = 100.;
akantu::Real dx = nodes(i,0) - length/4.;
akantu::Real dy = 0.0;
akantu::Real dz = 0.0;
if (spatial_dimension > 1) dy = nodes(i,1) - length/4.;
if (spatial_dimension == 3) dz = nodes(i,2) - length/4.;
akantu::Real d = sqrt(dx*dx + dy*dy + dz*dz);
// if(dx < 0.0){
if(d < 0.1){
boundary(i) = true;
temperature(i) = 300.;
}
}
-#ifdef AKANTU_USE_IOHELPER
iohelper::DumperParaview dumper;
paraviewInit(model,dumper);
-#endif
+
//main loop
int max_steps = 1000;
for(int i=0; i<max_steps; i++)
{
model.explicitPred();
model.updateResidual();
model.explicitCorr();
#ifdef AKANTU_USE_IOHELPER
if(i % 100 == 0)
paraviewDump(dumper);
#endif
if(i % 10 == 0)
std::cout << "Step " << i << "/" << max_steps << std::endl;
}
cout<< "\n\n Stable Time Step is : " << time_step << "\n \n" <<endl;
return 0;
}
-/* -------------------------------------------------------------------------- */
-#ifdef AKANTU_USE_IOHELPER
-/* -------------------------------------------------------------------------- */
void paraviewInit(akantu::HeatTransferModel & model, iohelper::Dumper & dumper) {
- // akantu::UInt nb_nodes = model.getFEM().getMesh().getNbNodes();
- // akantu::UInt nb_element = model.getFEM().getMesh().getNbElement(type);
+ base_name = "coordinates2";
+ akantu::UInt nb_nodes = model.getFEM().getMesh().getNbNodes();
+ akantu::UInt nb_element = model.getFEM().getMesh().getNbElement(type);
#pragma message "To change with new dumper"
- // // dumper.SetMode(iohelper::TEXT);
- // dumper.SetPoints(model.getFEM().getMesh().getNodes().values,
- // spatial_dimension, nb_nodes, "coordinates2");
- // dumper.SetConnectivity((int *)model.getFEM().getMesh().getConnectivity(type).values,
- // paraview_type, nb_element, iohelper::C_MODE);
- // dumper.AddNodeDataField(model.getTemperature().values,
- // 1, "temperature");
- // dumper.AddNodeDataField(model.getTemperatureRate().values,
+ // dumper.SetMode(iohelper::TEXT);
+ dumper.setPoints(model.getFEM().getMesh().getNodes().values,
+ spatial_dimension, nb_nodes, base_name);
+ dumper.setConnectivity((int *)model.getFEM().getMesh().getConnectivity(type).values,
+ paraview_type, nb_element, iohelper::C_MODE);
+ dumper.addNodeDataField("temperature",model.getTemperature().values,
+ 1,nb_nodes);
+ // dumper.addNodeDataField(model.getTemperatureRate().values,
// 1, "temperature_rate");
- // dumper.AddNodeDataField(model.getResidual().values,
+ // dumper.addNodeDataField(model.getResidual().values,
// 1, "residual");
- // dumper.AddNodeDataField(model.getCapacityLumped().values,
+ // dumper.addNodeDataField(model.getCapacityLumped().values,
// 1, "capacity_lumped");
- // dumper.AddElemDataField(model.getTemperatureGradient(type).values,
+ // dumper.addElemDataField(model.getTemperatureGradient(type).values,
// spatial_dimension, "temperature_gradient");
- // dumper.AddElemDataField(model.getConductivityOnQpoints(type).values,
+ // dumper.addElemDataField(model.getConductivityOnQpoints(type).values,
// spatial_dimension*spatial_dimension, "conductivity_qpoints");
- // dumper.SetPrefix("paraview/");
- // dumper.Init();
+ dumper.setPrefix("./");
+ dumper.init();
}
/* -------------------------------------------------------------------------- */
void paraviewDump(iohelper::Dumper & dumper) {
- // dumper.Dump();
-}
+ static int dump_step = 0;
+ std::stringstream sstr;
+ sstr << base_name << "-";
+ sstr.width(5);
+ sstr.fill('0');
+ sstr << dump_step;
-/* -------------------------------------------------------------------------- */
-#endif
-/* -------------------------------------------------------------------------- */
+ dumper.dump(sstr.str());
+
+ ++dump_step;
+}