diff --git a/src/io/parser/algebraic_parser.hh b/src/io/parser/algebraic_parser.hh
index b084f4696..cb0819e75 100644
--- a/src/io/parser/algebraic_parser.hh
+++ b/src/io/parser/algebraic_parser.hh
@@ -1,497 +1,497 @@
/**
* @file algebraic_parser.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Mon Jul 29 18:00:42 2013
*
* @brief algebraic_parser definition of the grammar
*
* @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/>.
*
*/
/* -------------------------------------------------------------------------- */
// Boost
#include <boost/config/warning_disable.hpp>
#include <boost/spirit/include/qi.hpp>
#include <boost/spirit/include/phoenix_core.hpp>
#include <boost/spirit/include/phoenix_operator.hpp>
#include <boost/spirit/include/phoenix_object.hpp>
#ifndef __AKANTU_ALGEBRAIC_PARSER_HH__
#define __AKANTU_ALGEBRAIC_PARSER_HH__
namespace spirit = boost::spirit;
namespace qi = boost::spirit::qi;
namespace lbs = boost::spirit::qi::labels;
namespace ascii = boost::spirit::ascii;
namespace phx = boost::phoenix;
namespace akantu {
namespace parser {
struct algebraic_error_handler_
{
template <typename, typename, typename>
struct result { typedef void type; };
template <typename Iterator>
void operator()(qi::info const& what,
Iterator err_pos,
Iterator last) const
{
AKANTU_EXCEPTION("Error! Expecting "
<< what // what failed?
<< " here: \""
<< std::string(err_pos, last) // iterators to error-pos, end
<< "\"");
}
};
static Real my_min(Real a, Real b) { return std::min(a, b); }
static Real my_max(Real a, Real b) { return std::max(a, b); }
template<class Iterator>
struct Skipper {
typedef qi::rule<Iterator, void()> type;
};
struct lazy_unary_func_ {
template <typename Funct, typename T>
struct result { typedef T type; };
template <typename Funct, typename T>
typename result<Funct, T>::type operator()(Funct f, T a) const {
return f(a);
}
};
struct lazy_binary_func_ {
template <typename Funct, typename T1, typename T2>
struct result { typedef T1 type; };
template <typename Funct, typename T1, typename T2>
typename result<Funct, T1, T2>::type operator()(Funct f, T1 a, T2 b) const {
return f(a, b);
}
};
struct lazy_pow_ {
template <typename T1, typename T2>
struct result { typedef T1 type; };
template <typename T1, typename T2>
typename result<T1, T2>::type operator()(T1 a, T2 b) const {
return std::pow(a, b);
}
};
struct lazy_eval_param_ {
template <typename T1, typename T2, typename res>
struct result { typedef res type; };
template <typename T1, typename T2, typename res>
res operator()(T1 a, const T2 & section,
__attribute__((unused)) const res & result) const {
return section.getParameter(a, _ppsc_current_and_parent_scope);
}
};
template<class Iterator, typename Skipper = spirit::unused_type>
struct AlgebraicGammar : qi::grammar<Iterator, Real(), Skipper> {
AlgebraicGammar(const ParserSection & section) : AlgebraicGammar::base_type(start,
"algebraic_grammar"),
section(section) {
phx::function<algebraic_error_handler_> const error_handler = algebraic_error_handler_();
phx::function<lazy_pow_> lazy_pow;
phx::function<lazy_unary_func_> lazy_unary_func;
phx::function<lazy_binary_func_> lazy_binary_func;
phx::function<lazy_eval_param_> lazy_eval_param;
start
= expr.alias()
;
expr
= term [ lbs::_val = lbs::_1 ]
>> *( ('+' > term [ lbs::_val += lbs::_1 ])
| ('-' > term [ lbs::_val -= lbs::_1 ])
)
;
term
= factor [ lbs::_val = lbs::_1 ]
>> *( ('*' > factor [ lbs::_val *= lbs::_1 ])
| ('/' > factor [ lbs::_val /= lbs::_1 ])
)
;
factor
= number [ lbs::_val = lbs::_1 ]
>> *("**" > number [ lbs::_val = lazy_pow(lbs::_val, lbs::_1) ])
;
number
= real [ lbs::_val = lbs::_1 ]
| ('-' > number [ lbs::_val = -lbs::_1 ])
| ('+' > number [ lbs::_val = lbs::_1 ])
| constant [ lbs::_val = lbs::_1 ]
| function [ lbs::_val = lbs::_1 ]
| ('(' > expr > ')') [ lbs::_val = lbs::_1 ]
| variable [ lbs::_val = lbs::_1 ]
;
function
= (qi::no_case[unary_function]
> '('
> expr
> ')') [ lbs::_val = lazy_unary_func(lbs::_1, lbs::_2) ]
| (qi::no_case[binary_function]
> '(' >> expr
> ',' >> expr
> ')') [ lbs::_val = lazy_binary_func(lbs::_1, lbs::_2, lbs::_3) ]
;
variable
= key [ lbs::_val = lazy_eval_param(lbs::_1, section, lbs::_val) ]
;
key
= qi::char_("a-zA-Z_") >> *qi::char_("a-zA-Z_0-9") // coming from the InputFileGrammar
;
constant.add
("pi", M_PI)
("e", M_E);
unary_function.add
("abs" , &std::abs )
("acos" , &std::acos )
("asin" , &std::asin )
("atan" , &std::atan )
("ceil" , &std::ceil )
("cos" , &std::cos )
("cosh" , &std::cosh )
("exp" , &std::exp )
("floor" , &std::floor )
("log10" , &std::log10 )
("log" , &std::log )
("sin" , &std::sin )
("sinh" , &std::sinh )
("sqrt" , &std::sqrt )
("tan" , &std::tan )
("tanh" , &std::tanh )
#if defined(AKANTU_CORE_CXX11 )
("acosh" , &std::acosh )
("asinh" , &std::asinh )
("atanh" , &std::atanh )
- ("crbt" , &std::crbt )
("exp2" , &std::exp2 )
("expm1" , &std::expm1 )
("log1p" , &std::log1p )
("log2" , &std::log2 )
("erf" , &std::erf )
("erfc" , &std::erfc )
("lgamma", &std::lgamma)
("tgamma", &std::tgamma)
("trunc" , &std::trunc )
("round" , &std::round )
+ // ("crbt" , &std::crbt )
#endif
;
binary_function.add
("pow" , &std::pow )
("min" , &parser::my_min)
("max" , &parser::my_max)
("atan2", &std::atan2 )
("fmod" , &std::fmod )
#if defined(AKANTU_CORE_CXX11)
("hypot", &std::hypot )
#endif
;
qi::on_error<qi::fail>(start, error_handler(lbs::_4, lbs::_3, lbs::_2));
expr .name("expression");
term .name("term");
factor .name("factor");
number .name("numerical-value");
variable.name("variable");
function.name("function");
constant.name("constants-list");
unary_function.name("unary-functions-list");
binary_function.name("binary-functions-list");
}
private:
qi::rule<Iterator, Real(), Skipper> start;
qi::rule<Iterator, Real(), Skipper> expr;
qi::rule<Iterator, Real(), Skipper> term;
qi::rule<Iterator, Real(), Skipper> factor;
qi::rule<Iterator, Real(), Skipper> number;
qi::rule<Iterator, Real(), Skipper> variable;
qi::rule<Iterator, Real(), Skipper> function;
qi::rule<Iterator, std::string(), Skipper> key;
qi::real_parser< Real, qi::real_policies<Real> > real;
qi::symbols<char, Real> constant;
qi::symbols<char, Real (*) (Real)> unary_function;
qi::symbols<char, Real (*) (Real, Real)> binary_function;
const ParserSection & section;
};
/* ---------------------------------------------------------------------- */
/* Vector Parser */
/* ---------------------------------------------------------------------- */
struct parsable_vector {
operator Vector<Real>() {
Vector<Real> tmp(_cells.size());
std::vector<Real>::iterator it = _cells.begin();
for (UInt i = 0; it != _cells.end(); ++it, ++i) tmp(i) = *it;
return tmp;
}
std::vector<Real> _cells;
};
struct parsable_matrix {
operator Matrix<Real>() {
size_t cols = 0;
std::vector<parsable_vector>::iterator it_rows = _cells.begin();
for (;it_rows != _cells.end(); ++it_rows) cols = std::max(cols, it_rows->_cells.size());
Matrix<Real> tmp(_cells.size(), _cells[0]._cells.size(), 0.);
it_rows = _cells.begin();
for (UInt i = 0; it_rows != _cells.end(); ++it_rows, ++i) {
std::vector<Real>::iterator it_cols = it_rows->_cells.begin();
for (UInt j = 0; it_cols != it_rows->_cells.end(); ++it_cols, ++j) {
tmp(i, j) = *it_cols;
}
}
return tmp;
}
std::vector<parsable_vector> _cells;
};
/* ---------------------------------------------------------------------- */
struct lazy_cont_add_ {
template <typename T1, typename T2>
struct result { typedef void type; };
template <typename T1, typename T2>
void operator()(T1 & cont, T2 & value) const {
cont._cells.push_back(value);
}
};
/* ---------------------------------------------------------------------- */
template<class Iterator, typename Skipper = spirit::unused_type>
struct VectorGammar : qi::grammar<Iterator, parsable_vector(), Skipper> {
VectorGammar(const ParserSection & section) : VectorGammar::base_type(start,
"vector_algebraic_grammar"),
number(section) {
phx::function<algebraic_error_handler_> const error_handler = algebraic_error_handler_();
//phx::function<lazy_algebraic_eval_> lazy_algebraic_eval;
phx::function<lazy_cont_add_> lazy_vector_add;
start
= '[' > vector > ']'
;
vector
= ( number [ lazy_vector_add(lbs::_a, lbs::_1) ]
>> *( ','
>> number [ lazy_vector_add(lbs::_a, lbs::_1) ]
)
) [ lbs::_val = lbs::_a ]
;
qi::on_error<qi::fail>(start, error_handler(lbs::_4, lbs::_3, lbs::_2));
start .name("start");
vector.name("vector");
number.name("value");
}
private:
qi::rule<Iterator, parsable_vector(), Skipper> start;
qi::rule<Iterator, parsable_vector(), qi::locals<parsable_vector>, Skipper> vector;
qi::rule<Iterator, Real(), Skipper> value;
AlgebraicGammar<Iterator, Skipper> number;
};
/* ---------------------------------------------------------------------- */
struct lazy_vector_eval_ {
template <typename T1, typename T2, typename res>
struct result { typedef bool type; };
template <typename T1, typename T2, typename res>
bool operator()(T1 a, const T2 & section, res & result) const {
std::string value = section.getParameter(a, _ppsc_current_and_parent_scope);
std::string::const_iterator b = value.begin();
std::string::const_iterator e = value.end();
parser::VectorGammar<std::string::const_iterator, qi::space_type> grammar(section);
return qi::phrase_parse(b, e, grammar, qi::space, result);
}
};
/* ---------------------------------------------------------------------- */
template<class Iterator, typename Skipper = spirit::unused_type>
struct MatrixGammar : qi::grammar<Iterator, parsable_matrix(), Skipper> {
MatrixGammar(const ParserSection & section) : MatrixGammar::base_type(start,
"matrix_algebraic_grammar"),
vector(section) {
phx::function<algebraic_error_handler_> const error_handler = algebraic_error_handler_();
phx::function<lazy_vector_eval_> lazy_vector_eval;
phx::function<lazy_cont_add_> lazy_matrix_add;
start
= '[' > matrix > ']'
;
matrix
= ( rows [ lazy_matrix_add(lbs::_a, lbs::_1) ]
>> *( ','
>> rows [ lazy_matrix_add(lbs::_a, lbs::_1) ]
)
) [ lbs::_val = lbs::_a ]
;
rows
= eval_vector
| vector
;
eval_vector
= (key [ lbs::_pass = lazy_vector_eval(lbs::_1, section, lbs::_a) ]) [lbs::_val = lbs::_a]
;
key
= qi::char_("a-zA-Z_") >> *qi::char_("a-zA-Z_0-9") // coming from the InputFileGrammar
;
qi::on_error<qi::fail>(start, error_handler(lbs::_4, lbs::_3, lbs::_2));
start.name("matrix");
matrix.name("all_rows");
rows .name("rows");
vector.name("vector");
eval_vector.name("eval_vector");
}
private:
qi::rule<Iterator, parsable_matrix(), Skipper> start;
qi::rule<Iterator, parsable_matrix(), qi::locals<parsable_matrix>, Skipper> matrix;
qi::rule<Iterator, parsable_vector(), Skipper> rows;
qi::rule<Iterator, parsable_vector(), qi::locals<parsable_vector>, Skipper> eval_vector;
qi::rule<Iterator, std::string(), Skipper> key;
VectorGammar<Iterator, Skipper> vector;
};
/* ---------------------------------------------------------------------- */
/* Randon Generator */
/* ---------------------------------------------------------------------- */
struct ParsableRandomGenerator {
ParsableRandomGenerator(Real base = Real(),
const RandomDistributionType & type = _rdt_not_defined,
const parsable_vector & parameters = parsable_vector()) :
base(base), type(type), parameters(parameters) {}
Real base;
RandomDistributionType type;
parsable_vector parameters;
};
/* ---------------------------------------------------------------------- */
template<class Iterator, typename Skipper = spirit::unused_type>
struct RandomGeneratorGammar : qi::grammar<Iterator, ParsableRandomGenerator(), Skipper> {
RandomGeneratorGammar(const ParserSection & section) : RandomGeneratorGammar::base_type(start,
"random_generator_grammar"),
number(section) {
phx::function<algebraic_error_handler_> const error_handler = algebraic_error_handler_();
phx::function<lazy_cont_add_> lazy_params_add;
start
= generator.alias()
;
generator
= qi::hold[distribution [ lbs::_val = lbs::_1 ] ]
| number [ lbs::_val = phx::construct<ParsableRandomGenerator>(lbs::_1) ]
;
distribution
= (
number
>> generator_type
>> '['
>> generator_params
>> ']'
) [ lbs::_val = phx::construct<ParsableRandomGenerator>(lbs::_1, lbs::_2, lbs::_3) ]
;
generator_params
= ( number [ lazy_params_add(lbs::_a, lbs::_1) ]
>> *( ','
> number [ lazy_params_add(lbs::_a, lbs::_1) ]
)
) [ lbs::_val = lbs::_a ]
;
#define AKANTU_RANDOM_DISTRIBUTION_TYPE_ADD(r, data, elem) \
(BOOST_PP_STRINGIZE(BOOST_PP_TUPLE_ELEM(2, 0, elem)), \
AKANTU_RANDOM_DISTRIBUTION_TYPES_PREFIX(BOOST_PP_TUPLE_ELEM(2, 0, elem)))
generator_type.add
BOOST_PP_SEQ_FOR_EACH(AKANTU_RANDOM_DISTRIBUTION_TYPE_ADD, _, AKANTU_RANDOM_DISTRIBUTION_TYPES);
#undef AKANTU_RANDOM_DISTRIBUTION_TYPE_ADD
qi::on_error<qi::fail>(start, error_handler(lbs::_4, lbs::_3, lbs::_2));
start .name("random-generator");
generator .name("random-generator");
distribution .name("random-distribution");
generator_type .name("generator-type");
generator_params.name("generator-parameters");
number .name("number");
}
private:
qi::rule<Iterator, ParsableRandomGenerator(), Skipper> start;
qi::rule<Iterator, ParsableRandomGenerator(), Skipper> generator;
qi::rule<Iterator, ParsableRandomGenerator(), Skipper> distribution;
qi::rule<Iterator, parsable_vector(), qi::locals<parsable_vector>, Skipper> generator_params;
AlgebraicGammar<Iterator, Skipper> number;
qi::symbols<char, RandomDistributionType> generator_type;
};
}
}
#endif /* __AKANTU_ALGEBRAIC_PARSER_HH__ */
diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh
index d58704a4c..97be2e14c 100644
--- a/src/model/solid_mechanics/solid_mechanics_model.hh
+++ b/src/model/solid_mechanics/solid_mechanics_model.hh
@@ -1,675 +1,675 @@
/**
* @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 "aka_types.hh"
#include "model.hh"
#include "data_accessor.hh"
#include "mesh.hh"
#include "dumpable.hh"
#include "boundary_condition.hh"
#include "integrator_gauss.hh"
#include "shape_lagrange.hh"
#include "integration_scheme_2nd_order.hh"
#include "solver.hh"
#include "material_selector.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
class Material;
class IntegrationScheme2ndOrder;
class Contact;
class SparseMatrix;
}
__BEGIN_AKANTU__
struct SolidMechanicsModelOptions : public ModelOptions {
SolidMechanicsModelOptions(AnalysisMethod analysis_method = _explicit_lumped_mass,
bool no_init_materials = false) :
analysis_method(analysis_method), no_init_materials(no_init_materials) {}
AnalysisMethod analysis_method;
bool no_init_materials;
};
extern const SolidMechanicsModelOptions default_solid_mechanics_model_options;
class SolidMechanicsModel : public Model, public DataAccessor,
public MeshEventHandler,
public Dumpable,
public BoundaryCondition<SolidMechanicsModel> {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
class NewMaterialElementsEvent : public NewElementsEvent {
public:
AKANTU_GET_MACRO_NOT_CONST(MaterialList, material, Array<UInt> &);
AKANTU_GET_MACRO(MaterialList, material, const Array<UInt> &);
protected:
Array<UInt> material;
};
typedef FEMTemplate<IntegratorGauss,ShapeLagrange> MyFEMType;
SolidMechanicsModel(Mesh & mesh,
UInt spatial_dimension = _all_dimensions,
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,
const ModelOptions & options = default_solid_mechanics_model_options);
/// initialize the fem object needed for boundary conditions
void initFEMBoundary();
/// register the tags associated with the parallel synchronizer
void initParallel(MeshPartition * partition, DataAccessor * data_accessor=NULL);
/// allocate all vectors
void initArrays();
/// allocate all vectors
void initArraysPreviousDisplacment();
/// initialize all internal arrays for materials
virtual 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(AnalysisMethod analysis_method = _explicit_lumped_mass);
bool isExplicit()
{ return method == _explicit_lumped_mass || method == _explicit_consistent_mass; }
/// 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
virtual 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();
void updateIncrement();
void updatePreviousDisplacement();
/// Solve the system @f[ A x = \alpha b @f] with A a lumped matrix
void solveLumped(Array<Real> & x,
const Array<Real> & A,
const Array<Real> & b,
const Array<bool> & boundary,
Real alpha);
/// explicit integration predictor
void explicitPred();
/// explicit integration corrector
void explicitCorr();
public:
void solveStep();
/* ------------------------------------------------------------------------ */
/* 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();
public:
template<SolveConvergenceMethod method, SolveConvergenceCriteria criteria>
bool solveStep(Real tolerance,
UInt max_iteration = 100);
public:
/// solve @f[ A\delta u = f_ext - f_int @f] in displacement
void solveDynamic();
/// solve Ku = f
void solveStatic();
/// solve Ku = f
void solveStatic(Array<bool> & boundary_normal, Array<Real> & EulerAngles);
/// test if the system is converged
template<SolveConvergenceCriteria criteria>
bool testConvergence(Real tolerance, Real & error);
/// test the convergence (norm of increment)
bool testConvergenceIncrement(Real tolerance) __attribute__((deprecated));
bool testConvergenceIncrement(Real tolerance, Real & error) __attribute__((deprecated));
/// test the convergence (norm of residual)
bool testConvergenceResidual(Real tolerance) __attribute__((deprecated));
bool testConvergenceResidual(Real tolerance, Real & error) __attribute__((deprecated));
/// create and return the velocity damping matrix
SparseMatrix & initVelocityDampingMatrix();
/// 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();
public:
//protected: Daniel changed it just for a test
/// compute A and solve @f[ A\delta u = f_ext - f_int @f]
template<NewmarkBeta::IntegrationSchemeCorrectorType type>
void solve(Array<Real> & increment);
/* ------------------------------------------------------------------------ */
/* Explicit/Implicit */
/* ------------------------------------------------------------------------ */
public:
/// Update the stresses for the computation of the residual of the Stiffness
/// matrix in the case of finite deformation
void computeStresses();
/// synchronize the ghost element boundaries values
void synchronizeBoundaries();
/* ------------------------------------------------------------------------ */
/* Materials (solid_mechanics_model_material.cc) */
/* ------------------------------------------------------------------------ */
public:
/// registers all the custom materials of a given type present in the input file
template <typename M>
void registerNewCustomMaterials(const ID & mat_type);
/// register an empty material of a given type
template <typename M>
Material & registerNewEmptyMaterial(const std::string & name);
// /// Use a UIntData in the mesh to specify the material to use per element
// void setMaterialIDsFromIntData(const std::string & data_name);
protected:
/// register a material in the dynamic database
template <typename M>
Material & registerNewMaterial(const ParserSection & mat_section);
/// read the material files to instantiate all the materials
void instantiateMaterials();
/* ------------------------------------------------------------------------ */
/* Mass (solid_mechanics_model_mass.cc) */
/* ------------------------------------------------------------------------ */
public:
/// assemble the lumped mass matrix
void assembleMassLumped();
/// assemble the mass matrix for consistent mass resolutions
void assembleMass();
protected:
/// assemble the lumped mass matrix for local and ghost elements
void assembleMassLumped(GhostType ghost_type);
/// assemble the mass matrix for either _ghost or _not_ghost elements
void assembleMass(GhostType ghost_type);
/// fill a vector of rho
void computeRho(Array<Real> & rho,
ElementType type,
GhostType ghost_type);
/* ------------------------------------------------------------------------ */
/* Data Accessor inherited members */
/* ------------------------------------------------------------------------ */
public:
inline virtual UInt getNbDataForElements(const Array<Element> & elements,
SynchronizationTag tag) const;
inline virtual void packElementData(CommunicationBuffer & buffer,
const Array<Element> & elements,
SynchronizationTag tag) const;
inline virtual void unpackElementData(CommunicationBuffer & buffer,
const Array<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 Array<Element> & elements,
Array<Element> * elements_per_mat) const;
inline virtual void packBarycenter(const Mesh & mesh,
CommunicationBuffer & buffer,
const Array<Element> & elements,
SynchronizationTag tag) const;
inline virtual void unpackBarycenter(const Mesh & mesh,
CommunicationBuffer & buffer,
const Array<Element> & elements,
SynchronizationTag tag);
/* ------------------------------------------------------------------------ */
/* Mesh Event Handler inherited members */
/* ------------------------------------------------------------------------ */
protected:
virtual void onNodesAdded (const Array<UInt> & nodes_list,
const NewNodesEvent & event);
virtual void onNodesRemoved(const Array<UInt> & element_list,
const Array<UInt> & new_numbering,
const RemovedNodesEvent & event);
virtual void onElementsAdded (const Array<Element> & nodes_list,
const NewElementsEvent & event);
virtual void onElementsRemoved(const Array<Element> & element_list,
const ByElementTypeUInt & new_numbering,
const RemovedElementsEvent & event);
/* ------------------------------------------------------------------------ */
/* Dumpable interface */
/* ------------------------------------------------------------------------ */
public:
virtual void addDumpFieldToDumper(const std::string & dumper_name,
const std::string & field_id);
virtual void addDumpGroupField(const std::string & field_id,
const std::string & group_name);
virtual void addDumpGroupFieldToDumper(const std::string & dumper_name,
const std::string & field_id,
const std::string & group_name);
virtual void removeDumpGroupField(const std::string & field_id,
const std::string & group_name);
virtual void removeDumpGroupFieldFromDumper(const std::string & dumper_name,
const std::string & field_id,
const std::string & group_name);
virtual void addDumpFieldVectorToDumper(const std::string & dumper_name,
const std::string & field_id);
virtual void addDumpGroupFieldVector(const std::string & field_id,
const std::string & group_name);
virtual void addDumpGroupFieldVectorToDumper(const std::string & dumper_name,
const std::string & field_id,
const std::string & group_name);
virtual void addDumpFieldTensorToDumper(const std::string & dumper_name,
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
void setTimeStep(Real time_step);
/// get the value of the conversion from forces/ mass to acceleration
AKANTU_GET_MACRO(F_M2A, f_m2a, Real);
/// set the value of the conversion from forces/ mass to acceleration
AKANTU_SET_MACRO(F_M2A, f_m2a, Real);
/// get the SolidMechanicsModel::displacement vector
AKANTU_GET_MACRO(Displacement, *displacement, Array<Real> &);
/// get the SolidMechanicsModel::previous_displacement vector
AKANTU_GET_MACRO(PreviousDisplacement, *previous_displacement, Array<Real> &);
/// get the SolidMechanicsModel::current_position vector \warn only consistent
/// after a call to SolidMechanicsModel::updateCurrentPosition
AKANTU_GET_MACRO(CurrentPosition, *current_position, const Array<Real> &);
/// get the SolidMechanicsModel::increment vector \warn only consistent if
/// SolidMechanicsModel::setIncrementFlagOn has been called before
AKANTU_GET_MACRO(Increment, *increment, Array<Real> &);
/// get the lumped SolidMechanicsModel::mass vector
AKANTU_GET_MACRO(Mass, *mass, Array<Real> &);
/// get the SolidMechanicsModel::velocity vector
AKANTU_GET_MACRO(Velocity, *velocity, Array<Real> &);
/// get the SolidMechanicsModel::acceleration vector, updated by
/// SolidMechanicsModel::updateAcceleration
AKANTU_GET_MACRO(Acceleration, *acceleration, Array<Real> &);
/// get the SolidMechanicsModel::force vector (boundary forces)
AKANTU_GET_MACRO(Force, *force, Array<Real> &);
/// get the SolidMechanicsModel::residual vector, computed by
/// SolidMechanicsModel::updateResidual
AKANTU_GET_MACRO(Residual, *residual, Array<Real> &);
/// get the SolidMechanicsModel::boundary vector
AKANTU_GET_MACRO(Boundary, *boundary, Array<bool> &);
/// get the SolidMechnicsModel::incrementAcceleration vector
AKANTU_GET_MACRO(IncrementAcceleration, *increment_acceleration, Array<Real> &);
/// get the value of the SolidMechanicsModel::increment_flag
AKANTU_GET_MACRO(IncrementFlag, increment_flag, bool);
/// get a particular material (by material index)
inline Material & getMaterial(UInt mat_index);
/// get a particular material (by material index)
inline const Material & getMaterial(UInt mat_index) const;
/// get a particular material (by material name)
inline Material & getMaterial(const std::string & name);
/// get a particular material (by material name)
inline const Material & getMaterial(const std::string & name) const;
/// get a particular material id from is name
inline UInt getMaterialIndex(const std::string & name) const;
/// give the number of materials
inline UInt getNbMaterials() const { return materials.size(); };
inline void setMaterialSelector(MaterialSelector & selector);
/// 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();
Real getKineticEnergy(const ElementType & type, UInt index);
/// 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, const 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 &);
/// get the velocity damping matrix
AKANTU_GET_MACRO(VelocityDampingMatrix, *velocity_damping_matrix, SparseMatrix &);
/// get the FEM object to integrate or interpolate on the boundary
inline FEM & getFEMBoundary(const ID & name = "");
/// get integrator
AKANTU_GET_MACRO(Integrator, *integrator, const IntegrationScheme2ndOrder &);
/// get access to the internal solver
AKANTU_GET_MACRO(Solver, *solver, Solver &);
/// get synchronizer
AKANTU_GET_MACRO(Synchronizer, *synch_parallel, const DistributedSynchronizer &);
AKANTU_GET_MACRO(ElementIndexByMaterial, element_index_by_material, const ByElementTypeArray<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);
protected:
- friend Material;
+ friend class Material;
template<UInt DIM, template <UInt> class WeightFunction>
friend class MaterialNonLocal;
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
Array<Real> * displacement;
/// displacements array at the previous time step (used in finite deformation)
Array<Real> * previous_displacement;
/// lumped mass array
Array<Real> * mass;
/// velocities array
Array<Real> * velocity;
/// accelerations array
Array<Real> * acceleration;
/// accelerations array
Array<Real> * increment_acceleration;
/// forces array
Array<Real> * force;
/// residuals array
Array<Real> * residual;
/// boundaries array
Array<bool> * boundary;
/// array of current position used during update residual
Array<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;
/// vectors containing local material element index for each global element index
ByElementTypeUInt element_index_by_material;
/// list of used materials
std::vector<Material *> materials;
/// mapping between material name and material internal id
std::map<std::string, UInt> materials_names_to_id;
/// class defining of to choose a material
MaterialSelector * material_selector;
/// define if it is the default selector or not
bool is_default_material_selector;
/// integration scheme of second order used
IntegrationScheme2ndOrder * integrator;
/// increment of displacement
Array<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;
/// analysis method check the list in akantu::AnalysisMethod
AnalysisMethod method;
/// internal synchronizer for parallel computations
DistributedSynchronizer * synch_parallel;
/// tells if the material are instantiated
bool are_materials_instantiated;
};
__END_AKANTU__
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#include "parser.hh"
#include "material.hh"
__BEGIN_AKANTU__
#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__
#include "material_selector_tmpl.hh"
#endif /* __AKANTU_SOLID_MECHANICS_MODEL_HH__ */