diff --git a/doc/manual/manual-io.tex b/doc/manual/manual-io.tex
index 51584d692..bc53574f1 100644
--- a/doc/manual/manual-io.tex
+++ b/doc/manual/manual-io.tex
@@ -1,105 +1,107 @@
\chapter{Input/Output}\index{I\/O}
\section{Generic data}
In this chapter, we address ways to get the internal data in human-readable formats.
The models in \akantu handle data associated to the
mesh, but this data can be split into several \code{Arrays}. For example, the
data stored per element type in a \code{ElementTypeMapArray} is composed of as
many \code{Array}s as types in the mesh.
In order to get this data in a visualization software, the models contain a
object to dump \code{VTK} files. These files can be visualized in software such
as \code{ParaView}\cite{paraview}, \code{ViSit}\cite{visit} or \code{Mayavi}\cite{mayavi}.
The internal dumper of the model can be configured to specify which data fields
are to be written. This is done with the
\code{addDumpField}\index{I\/O!addDumpField} method. By default all the files
are generated in a folder called \code{paraview/}
\begin{cpp}
model.setBaseName("output"); // prefix for all generated files
model.addDumpField("displacement");
model.addDumpField("stress");
...
model.dump()
\end{cpp}
The fields are dumped with the number of components of the memory. For example, in 2D, the memory has
\code{Vector}s of 2 components, or the $2^{nd}$ order tensors with $2\times2$ components.
This memory can be dealt with \code{addDumpFieldVector}\index{I\/O!addDumpFieldVector} which always dumps
\code{Vector}s with 3 components or \code{addDumpFieldTensor}\index{I\/O!addDumpFieldTensor} which dumps $2^{nd}$
order tensors with $3\times3$ components respectively. The routines \code{addDumpFieldVector}\index{I\/O!addDumpFieldVector} and
\code{addDumpFieldTensor}\index{I\/O!addDumpFieldTensor} were introduced because of Paraview which mostly manipulate 3D
data.
Those fields which are stored by quadrature point are modified to be seen in the
\code{VTK} file as elemental data. To do this, the default is to average the
values of all the quadrature points.
The list of fields depends on the models (for
\code{SolidMechanicsModel} see table~\ref{tab:io:smm_field_list}).
\begin{table}
\centering
\begin{tabular}{llll}
\toprule
key & type & support \\
\midrule
displacement & Vector<Real> & nodes \\
velocity & Vector<Real> & nodes \\
acceleration & Vector<Real> & nodes \\
force & Vector<Real> & nodes \\
residual & Vector<Real> & nodes \\
boundary & Vector<bool> & nodes \\
mass & Vector<Real> & nodes \\
partitions & Real & elements \\
stress & Matrix<Real> & quadrature points \\
+ Von Mises stress & Real & quadrature points \\
+ grad_u & Matrix<Real> & quadrature points \\
strain & Matrix<Real> & quadrature points \\
\textit{material internals} & variable & quadrature points \\
\bottomrule
\end{tabular}
\caption{List of dumpable fields for \code{SolidMechanicsModel}.}
\label{tab:io:smm_field_list}
\end{table}
The user can also register external fields which have the same mesh as the mesh from the model as support. To do this, an object of type \code{Field} has to be created.\index{I\/O!addDumpFieldExternal}
\begin{itemize}
\item For nodal fields :
\begin{cpp}
Vector<T> vect(nb_nodes, nb_component);
Field field = new DumperIOHelper::NodalField<T>(vect));
model.addDumpFieldExternal("my_field", field);
\end{cpp}
\item For elemental fields :
\begin{cpp}
ElementTypeMapArray<T> arr;
Field field = new DumperIOHelper::ElementalField<T>(arr, spatial_displacement));
model.addDumpFieldExternal("my_field", field);
\end{cpp}
\end{itemize}
\section{Cohesive elements' data}
Cohesive elements and their relative data can be easily dumped thanks
to a specific dumper contained in
\code{SolidMechanicsModelCohesive}. In order to use it, one has just
to add the string \code{"cohesive elements"} when calling each method
already illustrated. Here is an example on how to dump displacement
and damage:
\begin{cpp}
model.setBaseNameToDumper("cohesive elements", "cohesive_elements_output");
model.addDumpFieldVectorToDumper("cohesive elements", "displacement");
model.addDumpFieldToDumper("cohesive elements", "damage");
...
model.dump("cohesive elements");
\end{cpp}
%%% Local Variables:
%%% mode: latex
%%% TeX-master: "manual"
%%% End:
diff --git a/src/io/dumper/dumper_material_padders.hh b/src/io/dumper/dumper_material_padders.hh
index db4a43bb2..67eae531d 100644
--- a/src/io/dumper/dumper_material_padders.hh
+++ b/src/io/dumper/dumper_material_padders.hh
@@ -1,167 +1,237 @@
/**
* @file dumper_material_padders.hh
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Tue Sep 02 2014
* @date last modification: Fri Sep 19 2014
*
* @brief Material padders for plane stress/ plane strain
*
* @section LICENSE
*
* Copyright (©) 2014 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_DUMPER_MATERIAL_PADDERS_HH__
#define __AKANTU_DUMPER_MATERIAL_PADDERS_HH__
/* -------------------------------------------------------------------------- */
#include "dumper_padding_helper.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
__BEGIN_AKANTU_DUMPER__
/* -------------------------------------------------------------------------- */
-
-
-template<class T, class R>
-class MaterialPadder : public PadderGeneric<Vector<T>, R > {
-
+class MaterialFunctor {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
-
public:
- MaterialPadder(const SolidMechanicsModel & model) :
+ MaterialFunctor(const SolidMechanicsModel & model) :
model(model),
- material_index(model.getMaterialByElement()) {}
+ material_index(model.getMaterialByElement()),
+ spatial_dimension(model.getSpatialDimension()){}
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
-
/// return the material from the global element index
const Material & getMaterialFromGlobalIndex(Element global_index){
-
UInt index = global_index.getIndex();
UInt material_id = material_index(global_index.getType())(index);
const Material & material = model.getMaterial(material_id);
return material;
}
/// return the type of the element from global index
ElementType getElementTypeFromGlobalIndex(Element global_index){
return global_index.getType();
}
protected:
-
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
/// all material padders probably need access to solid mechanics model
const SolidMechanicsModel & model;
+
/// they also need an access to the map from global ids to material id and local ids
const ElementTypeMapArray<UInt> & material_index;
/// the number of data per element
const ElementTypeMapArray<UInt> nb_data_per_element;
+ UInt spatial_dimension;
+};
+
+/* -------------------------------------------------------------------------- */
+template<class T, class R>
+class MaterialPadder : public MaterialFunctor, public PadderGeneric<Vector<T>, R > {
+public:
+ MaterialPadder(const SolidMechanicsModel & model) :
+ MaterialFunctor(model) {}
};
/* -------------------------------------------------------------------------- */
template <UInt spatial_dimension>
-class StressPadder :
- public MaterialPadder<Real,Matrix<Real> > {
+class StressPadder : public MaterialPadder<Real, Matrix<Real> > {
public:
StressPadder(const SolidMechanicsModel & model) :
MaterialPadder<Real, Matrix<Real> >(model){
this->setPadding(3,3);
}
inline Matrix<Real> func(const Vector<Real> & in, Element global_element_id){
UInt nrows = spatial_dimension;
UInt ncols = in.size() / nrows;
UInt nb_data = in.size() / (nrows*nrows);
- Matrix<Real> stress = this->pad(in, nrows,ncols, nb_data);
+ Matrix<Real> stress = this->pad(in, nrows, ncols, nb_data);
const Material & material = this->getMaterialFromGlobalIndex(global_element_id);
bool plane_strain = true;
if(spatial_dimension == 2)
plane_strain = !material.getParam<bool>("Plane_Stress");
if(plane_strain) {
Real nu = material.getParam<Real>("nu");
for (UInt d = 0; d < nb_data; ++d) {
stress(2, 2 + 3*d) = nu * (stress(0, 0 + 3*d) + stress(1, 1 + 3*d));
}
}
return stress;
}
UInt getDim(){return 9;};
UInt getNbComponent(UInt old_nb_comp){
return this->getDim();
};
};
/* -------------------------------------------------------------------------- */
template<UInt spatial_dimension>
-class StrainPadder : public MaterialPadder<Real, Matrix<Real> > {
+class StrainPadder : public MaterialFunctor, public PadderGeneric< Matrix<Real>, Matrix<Real> > {
public:
StrainPadder(const SolidMechanicsModel & model) :
- MaterialPadder<Real, Matrix<Real> >(model) {
+ MaterialFunctor(model) {
this->setPadding(3,3);
}
- inline Matrix<Real> func(const Vector<Real> & in, Element global_element_id){
-
+ inline Matrix<Real> func(const Matrix<Real> & in, Element global_element_id){
UInt nrows = spatial_dimension;
- UInt ncols = in.size() / nrows;
- UInt nb_data = in.size() / ncols;
+ UInt nb_data = in.size() / (nrows*nrows);
- Matrix<Real> strain = this->pad(in, nrows,ncols, nb_data);
+ Matrix<Real> strain = this->pad(in, nb_data);
const Material & material = this->getMaterialFromGlobalIndex(global_element_id);
bool plane_stress = material.getParam<bool>("Plane_Stress");
+
if(plane_stress) {
Real nu = material.getParam<Real>("nu");
for (UInt d = 0; d < nb_data; ++d) {
strain(2, 2 + 3*d) = nu / (nu - 1) * (strain(0, 0 + 3*d) + strain(1, 1 + 3*d));
}
}
+
return strain;
}
UInt getDim(){return 9;};
UInt getNbComponent(UInt old_nb_comp){
return this->getDim();
};
};
+/* -------------------------------------------------------------------------- */
+class ComputeStrain : public MaterialFunctor, public ComputeFunctor<Vector<Real>, Matrix<Real> > {
+public:
+ ComputeStrain(const SolidMechanicsModel & model) : MaterialFunctor(model) { }
+
+ inline Matrix<Real> func(const Vector<Real> & in, Element global_element_id){
+ UInt nrows = spatial_dimension;
+ UInt ncols = in.size() / nrows;
+ UInt nb_data = in.size() / (nrows*nrows);
+
+ Matrix<Real> ret_all_strain(nrows, ncols);
+ Tensor3<Real> all_grad_u(in.storage(), nrows, nrows, nb_data);
+ Tensor3<Real> all_strain(ret_all_strain.storage(), nrows, nrows, nb_data);
+
+ for (UInt d = 0; d < nb_data; ++d) {
+ Matrix<Real> grad_u = all_grad_u(d);
+ Matrix<Real> strain = all_strain(d);
+
+ if (spatial_dimension == 2)
+ Material::gradUToEpsilon<2>(grad_u, strain);
+ else if (spatial_dimension == 3)
+ Material::gradUToEpsilon<3>(grad_u, strain);
+ }
+
+ return ret_all_strain;
+ }
+
+ UInt getDim() { return spatial_dimension*spatial_dimension; };
+
+ UInt getNbComponent(UInt old_nb_comp){
+ return this->getDim();
+ };
+
+};
+
+/* -------------------------------------------------------------------------- */
+
+/* -------------------------------------------------------------------------- */
+class ComputeVonMisesStress : public MaterialFunctor,
+ public ComputeFunctor<Vector<Real>, Vector<Real> > {
+public:
+ ComputeVonMisesStress(const SolidMechanicsModel & model) : MaterialFunctor(model) { }
+
+ inline Vector<Real> func(const Vector<Real> & in, Element global_element_id){
+ UInt nrows = spatial_dimension;
+ UInt nb_data = in.size() / (nrows*nrows);
+
+ Vector<Real> von_mises_stress(nb_data);
+ Matrix<Real> deviatoric_stress(3, 3);
+
+ for (UInt d = 0; d < nb_data; ++d) {
+ Matrix<Real> cauchy_stress(in.storage() + d * nrows*nrows, nrows, nrows);
+ von_mises_stress(d) = Material::stressToVonMises(cauchy_stress);
+ }
+
+ return von_mises_stress;
+ }
+
+ UInt getDim() { return 1; };
+
+ UInt getNbComponent(UInt old_nb_comp){
+ return this->getDim();
+ };
+
+};
+
+
+
/* -------------------------------------------------------------------------- */
__END_AKANTU_DUMPER__
__END_AKANTU__
#endif /* __AKANTU_DUMPER_MATERIAL_PADDERS_HH__ */
diff --git a/src/model/solid_mechanics/material.hh b/src/model/solid_mechanics/material.hh
index b53e92ae8..ac55f055e 100644
--- a/src/model/solid_mechanics/material.hh
+++ b/src/model/solid_mechanics/material.hh
@@ -1,596 +1,600 @@
/**
* @file material.hh
*
* @author Marco Vocialta <marco.vocialta@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
*
* @date creation: Tue Jul 27 2010
* @date last modification: Tue Sep 16 2014
*
* @brief Mother class for all materials
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "aka_memory.hh"
#include "aka_voigthelper.hh"
#include "parser.hh"
#include "parsable.hh"
#include "data_accessor.hh"
#include "internal_field.hh"
#include "random_internal_field.hh"
#include "solid_mechanics_model_event_handler.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MATERIAL_HH__
#define __AKANTU_MATERIAL_HH__
/* -------------------------------------------------------------------------- */
namespace akantu {
class Model;
class SolidMechanicsModel;
}
__BEGIN_AKANTU__
/**
* Interface of all materials
* Prerequisites for a new material
* - inherit from this class
* - implement the following methods:
* \code
* virtual Real getStableTimeStep(Real h, const Element & element = ElementNull);
*
* virtual void computeStress(ElementType el_type,
* GhostType ghost_type = _not_ghost);
*
* virtual void computeTangentStiffness(const ElementType & el_type,
* Array<Real> & tangent_matrix,
* GhostType ghost_type = _not_ghost);
* \endcode
*
*/
class Material : public Memory, public DataAccessor, public Parsable,
public MeshEventHandler,
protected SolidMechanicsModelEventHandler {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
Material(SolidMechanicsModel & model, const ID & id = "");
virtual ~Material();
/* ------------------------------------------------------------------------ */
/* Function that materials can/should reimplement */
/* ------------------------------------------------------------------------ */
protected:
/// constitutive law
virtual void computeStress(__attribute__((unused)) ElementType el_type,
__attribute__((unused)) GhostType ghost_type = _not_ghost) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/// compute the tangent stiffness matrix
virtual void computeTangentModuli(__attribute__((unused)) const ElementType & el_type,
__attribute__((unused)) Array<Real> & tangent_matrix,
__attribute__((unused)) GhostType ghost_type = _not_ghost) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/// compute the potential energy
virtual void computePotentialEnergy(ElementType el_type, GhostType ghost_type = _not_ghost);
/// compute the potential energy for an element
virtual void computePotentialEnergyByElement(__attribute__((unused)) ElementType type,
__attribute__((unused)) UInt index,
__attribute__((unused)) Vector<Real> & epot_on_quad_points) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
virtual void updateEnergies(__attribute__((unused)) ElementType el_type,
__attribute__((unused)) GhostType ghost_type = _not_ghost) { }
virtual void updateEnergiesAfterDamage(__attribute__((unused)) ElementType el_type,
__attribute__((unused)) GhostType ghost_type = _not_ghost) {}
/// set the material to steady state (to be implemented for materials that need it)
virtual void setToSteadyState(__attribute__((unused)) ElementType el_type,
__attribute__((unused)) GhostType ghost_type = _not_ghost) { }
/// function called to update the internal parameters when the modifiable
/// parameters are modified
virtual void updateInternalParameters() {}
public:
/// compute the p-wave speed in the material
virtual Real getPushWaveSpeed(const Element & element) const { AKANTU_DEBUG_TO_IMPLEMENT(); }
/// compute the s-wave speed in the material
virtual Real getShearWaveSpeed(const Element & element) const { AKANTU_DEBUG_TO_IMPLEMENT(); }
/// get a material celerity to compute the stable time step (default: is the push wave speed)
virtual Real getCelerity(const Element & element) const { return getPushWaveSpeed(element); }
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
template<typename T>
void registerInternal(__attribute__((unused)) InternalField<T> & vect) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
template<typename T>
void unregisterInternal(__attribute__((unused)) InternalField<T> & vect) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/// initialize the material computed parameter
virtual void initMaterial();
/// compute the residual for this material
virtual void updateResidual(GhostType ghost_type = _not_ghost);
/// assemble the residual for this material
virtual void assembleResidual(GhostType ghost_type);
/// Operations before and after solveStep in implicit
virtual void beforeSolveStep() {}
virtual void afterSolveStep() {}
/// save the stress in the previous_stress if needed
virtual void savePreviousState();
/// compute the stresses for this material
virtual void computeAllStresses(GhostType ghost_type = _not_ghost);
virtual void computeAllNonLocalStresses(__attribute__((unused)) GhostType ghost_type = _not_ghost) {};
virtual void computeAllStressesFromTangentModuli(GhostType ghost_type = _not_ghost);
virtual void computeAllCauchyStresses(GhostType ghost_type = _not_ghost);
/// set material to steady state
void setToSteadyState(GhostType ghost_type = _not_ghost);
/// compute the stiffness matrix
virtual void assembleStiffnessMatrix(GhostType ghost_type);
/// add an element to the local mesh filter
inline UInt addElement(const ElementType & type,
UInt element,
const GhostType & ghost_type);
/// add many elements at once
void addElements(const Array<Element> & elements_to_add);
/// remove many element at once
void removeElements(const Array<Element> & elements_to_remove);
/// function to print the contain of the class
virtual void printself(std::ostream & stream, int indent = 0) const;
/**
* interpolate stress on given positions for each element by means
* of a geometrical interpolation on quadrature points
*/
void interpolateStress(ElementTypeMapArray<Real> & result,
const GhostType ghost_type = _not_ghost);
/**
* interpolate stress on given positions for each element by means
* of a geometrical interpolation on quadrature points and store the
* results per facet
*/
void interpolateStressOnFacets(ElementTypeMapArray<Real> & result,
const GhostType ghost_type = _not_ghost);
/**
* function to initialize the elemental field interpolation
* function by inverting the quadrature points' coordinates
*/
void initElementalFieldInterpolation(const ElementTypeMapArray<Real> & interpolation_points_coordinates);
/* ------------------------------------------------------------------------ */
/* Common part */
/* ------------------------------------------------------------------------ */
protected:
/// assemble the residual
template<UInt dim>
void assembleResidual(GhostType ghost_type);
/// Computation of Cauchy stress tensor in the case of finite deformation
template<UInt dim>
void computeCauchyStress(__attribute__((unused)) ElementType el_type,
__attribute__((unused)) GhostType ghost_type = _not_ghost);
template<UInt dim >
inline void computeCauchyStressOnQuad(const Matrix<Real> & F, const Matrix<Real> & S,
Matrix<Real> & cauchy,
const Real & C33 = 1.0 ) const;
template<UInt dim>
void computeAllStressesFromTangentModuli(const ElementType & type,
GhostType ghost_type);
template<UInt dim>
void assembleStiffnessMatrix(const ElementType & type,
GhostType ghost_type);
/// assembling in finite deformation
template<UInt dim>
void assembleStiffnessMatrixNL(const ElementType & type,
GhostType ghost_type);
template<UInt dim>
void assembleStiffnessMatrixL2(const ElementType & type,
GhostType ghost_type);
/// write the stress tensor in the Voigt notation.
template<UInt dim>
inline void SetCauchyStressArray(const Matrix<Real> & S_t, Matrix<Real> & Stress_vect);
inline UInt getTangentStiffnessVoigtSize(UInt spatial_dimension) const;
/// Size of the Stress matrix for the case of finite deformation see: Bathe et al, IJNME, Vol 9, 353-386, 1975
inline UInt getCauchyStressMatrixSize(UInt spatial_dimension) const;
/// Sets the stress matrix according to Bathe et al, IJNME, Vol 9, 353-386, 1975
template<UInt dim>
inline void setCauchyStressMatrix(const Matrix<Real> & S_t,
Matrix<Real> & Stress_matrix);
/// compute the potential energy by element
void computePotentialEnergyByElements();
/// resize the intenals arrays
void resizeInternals();
public:
/// compute the coordinates of the quadrature points
void computeQuadraturePointsCoordinates(ElementTypeMapArray<Real> & quadrature_points_coordinates,
const GhostType & ghost_type) const;
protected:
/// interpolate an elemental field on given points for each element
void interpolateElementalField(const ElementTypeMapArray<Real> & field,
ElementTypeMapArray<Real> & result,
const GhostType ghost_type);
/// interpolate an elemental field and store the results per facet
void interpolateElementalFieldOnFacets(const ElementTypeMapArray<Real> & field,
ElementTypeMapArray<Real> & result,
const GhostType ghost_type);
/// template function to initialize the elemental field interpolation
template <ElementType type>
void initElementalFieldInterpolation(const Array<Real> & quad_coordinates,
const Array<Real> & interpolation_points_coordinates,
const UInt nb_interpolation_points_per_elem,
const GhostType ghost_type);
/// build the coordinate matrix for the interpolation on elemental field
template <ElementType type>
inline void buildElementalFieldInterpolationCoodinates(const Matrix<Real> & coordinates,
Matrix<Real> & coordMatrix);
/// build interpolation coordinates for basic linear elements
inline void buildElementalFieldInterpolationCoodinatesLinear(const Matrix<Real> & coordinates,
Matrix<Real> & coordMatrix);
/// build interpolation coordinates for basic quadratic elements
inline void buildElementalFieldInterpolationCoodinatesQuadratic(const Matrix<Real> & coordinates,
Matrix<Real> & coordMatrix);
public:
/* ------------------------------------------------------------------------ */
/* Conversion functions */
/* ------------------------------------------------------------------------ */
template<UInt dim>
- inline void gradUToF (const Matrix<Real> & grad_u, Matrix<Real> & F) const;
- inline void rightCauchy(const Matrix<Real> & F, Matrix<Real> & C) const;
- inline void leftCauchy (const Matrix<Real> & F, Matrix<Real> & B) const;
+ static inline void gradUToF (const Matrix<Real> & grad_u, Matrix<Real> & F);
+ static inline void rightCauchy(const Matrix<Real> & F, Matrix<Real> & C);
+ static inline void leftCauchy (const Matrix<Real> & F, Matrix<Real> & B);
template<UInt dim>
- inline void gradUToEpsilon(const Matrix<Real> & grad_u, Matrix<Real> & epsilon) const;
+ static inline void gradUToEpsilon(const Matrix<Real> & grad_u,
+ Matrix<Real> & epsilon);
template<UInt dim>
- inline void gradUToGreenStrain(const Matrix<Real> & grad_u,
- Matrix<Real> & epsilon) const;
+ static inline void gradUToGreenStrain(const Matrix<Real> & grad_u,
+ Matrix<Real> & epsilon);
+
+ static inline Real stressToVonMises(const Matrix<Real> & stress);
+
protected:
/// converts global element to local element
inline Element convertToLocalElement(const Element & global_element) const;
/// converts local element to global element
inline Element convertToGlobalElement(const Element & local_element) const;
/// converts global quadrature point to local quadrature point
inline QuadraturePoint convertToLocalPoint(const QuadraturePoint & global_point) const;
/// converts local quadrature point to global quadrature point
inline QuadraturePoint convertToGlobalPoint(const QuadraturePoint & local_point) const;
/* ------------------------------------------------------------------------ */
/* DataAccessor inherited members */
/* ------------------------------------------------------------------------ */
public:
virtual inline UInt getNbDataForElements(const Array<Element> & elements,
SynchronizationTag tag) const;
virtual inline void packElementData(CommunicationBuffer & buffer,
const Array<Element> & elements,
SynchronizationTag tag) const;
virtual inline void unpackElementData(CommunicationBuffer & buffer,
const Array<Element> & elements,
SynchronizationTag tag);
template<typename T>
inline void packElementDataHelper(const ElementTypeMapArray<T> & data_to_pack,
CommunicationBuffer & buffer,
const Array<Element> & elements,
const ID & fem_id = ID()) const;
template<typename T>
inline void unpackElementDataHelper(ElementTypeMapArray<T> & data_to_unpack,
CommunicationBuffer & buffer,
const Array<Element> & elements,
const ID & fem_id = ID());
/* ------------------------------------------------------------------------ */
/* MeshEventHandler inherited members */
/* ------------------------------------------------------------------------ */
public:
/* ------------------------------------------------------------------------ */
virtual void onElementsAdded(const Array<Element> & element_list,
const NewElementsEvent & event);
virtual void onElementsRemoved(const Array<Element> & element_list,
const ElementTypeMapArray<UInt> & new_numbering,
const RemovedElementsEvent & event);
/* ------------------------------------------------------------------------ */
/* SolidMechanicsModelEventHandler inherited members */
/* ------------------------------------------------------------------------ */
public:
virtual void onBeginningSolveStep(const AnalysisMethod & method);
virtual void onEndSolveStep(const AnalysisMethod & method);
virtual void onDamageIteration();
virtual void onDamageUpdate();
virtual void onDump();
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
AKANTU_GET_MACRO(Model, *model, const SolidMechanicsModel &)
AKANTU_GET_MACRO(ID, Memory::getID(), const ID &);
AKANTU_GET_MACRO(Rho, rho, Real);
AKANTU_SET_MACRO(Rho, rho, Real);
/// return the potential energy for the subset of elements contained by the material
Real getPotentialEnergy();
/// return the potential energy for the provided element
Real getPotentialEnergy(ElementType & type, UInt index);
/// return the energy (identified by id) for the subset of elements contained by the material
virtual Real getEnergy(std::string energy_id);
/// return the energy (identified by id) for the provided element
virtual Real getEnergy(std::string energy_id, ElementType type, UInt index);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementFilter, element_filter, UInt);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(GradU, gradu, Real);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real);
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(PotentialEnergy, potential_energy, Real);
AKANTU_GET_MACRO(GradU, gradu, const ElementTypeMapArray<Real> &);
AKANTU_GET_MACRO(Stress, stress, const ElementTypeMapArray<Real> &);
AKANTU_GET_MACRO(ElementFilter, element_filter, const ElementTypeMapArray<UInt> &);
bool isNonLocal() const { return is_non_local; }
const Array<Real> & getArray(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost) const;
Array<Real> & getArray(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost);
const InternalField<Real> & getInternal(const ID & id) const;
InternalField<Real> & getInternal(const ID & id);
inline bool isInternal(const ID & id, const ElementKind & element_kind) const;
inline ElementTypeMap<UInt> getInternalDataPerElem(const ID & id, const ElementKind & element_kind) const;
bool isFiniteDeformation() const { return finite_deformation; }
bool isInelasticDeformation() const { return inelastic_deformation; }
template <typename T>
inline void setParam(const ID & param, T value);
template <typename T>
inline const T & getParam(const ID & param) const;
void flattenInternal(const std::string & field_id,
ElementTypeMapArray<Real> & internal_flat,
const GhostType ghost_type = _not_ghost,
ElementKind element_kind = _ek_not_defined);
protected:
bool isInit() const { return is_init; }
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
/// boolean to know if the material has been initialized
bool is_init;
std::map<ID, InternalField<Real> *> internal_vectors_real;
std::map<ID, InternalField<UInt> *> internal_vectors_uint;
protected:
/// Finite deformation
bool finite_deformation;
/// Finite deformation
bool inelastic_deformation;
/// material name
std::string name;
/// The model to witch the material belong
SolidMechanicsModel * model;
/// density : rho
Real rho;
/// spatial dimension
UInt spatial_dimension;
/// list of element handled by the material
ElementTypeMapArray<UInt> element_filter;
/// stresses arrays ordered by element types
InternalField<Real> stress;
/// eigenstrain arrays ordered by element types
InternalField<Real> eigenstrain;
/// grad_u arrays ordered by element types
InternalField<Real> gradu;
/// Green Lagrange strain (Finite deformation)
InternalField<Real> green_strain;
/// Second Piola-Kirchhoff stress tensor arrays ordered by element types (Finite deformation)
InternalField<Real> piola_kirchhoff_2;
/// potential energy by element
InternalField<Real> potential_energy;
/// tell if using in non local mode or not
bool is_non_local;
/// tell if the material need the previous stress state
bool use_previous_stress;
/// tell if the material need the previous strain state
bool use_previous_gradu;
/// elemental field interpolation coordinates
InternalField<Real> interpolation_inverse_coordinates;
/// elemental field interpolation points
InternalField<Real> interpolation_points_matrices;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#include "material_inline_impl.cc"
/// standard output stream operator
inline std::ostream & operator <<(std::ostream & stream, const Material & _this)
{
_this.printself(stream);
return stream;
}
__END_AKANTU__
#include "internal_field_tmpl.hh"
#include "random_internal_field_tmpl.hh"
/* -------------------------------------------------------------------------- */
/* Auto loop */
/* -------------------------------------------------------------------------- */
#define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type) \
Array<Real>::matrix_iterator gradu_it = \
this->gradu(el_type, ghost_type).begin(this->spatial_dimension, \
this->spatial_dimension); \
Array<Real>::matrix_iterator gradu_end = \
this->gradu(el_type, ghost_type).end(this->spatial_dimension, \
this->spatial_dimension); \
\
this->stress(el_type, \
ghost_type).resize(this->gradu(el_type, \
ghost_type).getSize()); \
\
Array<Real>::iterator< Matrix<Real> > stress_it = \
this->stress(el_type, ghost_type).begin(this->spatial_dimension, \
this->spatial_dimension); \
\
if(this->isFiniteDeformation()){ \
this->piola_kirchhoff_2(el_type, \
ghost_type).resize(this->gradu(el_type, \
ghost_type).getSize()); \
stress_it = \
this->piola_kirchhoff_2(el_type, \
ghost_type).begin(this->spatial_dimension, \
this->spatial_dimension); \
} \
\
for(;gradu_it != gradu_end; ++gradu_it, ++stress_it) { \
Matrix<Real> & __attribute__((unused)) grad_u = *gradu_it; \
Matrix<Real> & __attribute__((unused)) sigma = *stress_it
#define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END \
} \
#define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_mat) \
Array<Real>::matrix_iterator gradu_it = \
this->gradu(el_type, ghost_type).begin(this->spatial_dimension, \
this->spatial_dimension); \
Array<Real>::matrix_iterator gradu_end = \
this->gradu(el_type, ghost_type).end(this->spatial_dimension, \
this->spatial_dimension); \
Array<Real>::matrix_iterator sigma_it = \
this->stress(el_type, ghost_type).begin(this->spatial_dimension, \
this->spatial_dimension); \
\
tangent_mat.resize(this->gradu(el_type, ghost_type).getSize()); \
\
UInt tangent_size = \
this->getTangentStiffnessVoigtSize(this->spatial_dimension); \
Array<Real>::matrix_iterator tangent_it = \
tangent_mat.begin(tangent_size, \
tangent_size); \
\
for(;gradu_it != gradu_end; ++gradu_it, ++sigma_it, ++tangent_it) { \
Matrix<Real> & __attribute__((unused)) grad_u = *gradu_it; \
Matrix<Real> & __attribute__((unused)) sigma_tensor = *sigma_it; \
Matrix<Real> & tangent = *tangent_it
#define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END \
} \
/* -------------------------------------------------------------------------- */
#define INSTANSIATE_MATERIAL(mat_name) \
template class mat_name<1>; \
template class mat_name<2>; \
template class mat_name<3>
#endif /* __AKANTU_MATERIAL_HH__ */
diff --git a/src/model/solid_mechanics/material_inline_impl.cc b/src/model/solid_mechanics/material_inline_impl.cc
index 9819833f6..28473d8a4 100644
--- a/src/model/solid_mechanics/material_inline_impl.cc
+++ b/src/model/solid_mechanics/material_inline_impl.cc
@@ -1,476 +1,490 @@
/**
* @file material_inline_impl.cc
*
* @author Marco Vocialta <marco.vocialta@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
*
* @date creation: Tue Jul 27 2010
* @date last modification: Tue Sep 16 2014
*
* @brief Implementation of the inline functions of the class material
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014 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/>.
*
*/
__END_AKANTU__
#include "solid_mechanics_model.hh"
#include <iostream>
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
inline UInt Material::addElement(const ElementType & type,
UInt element,
const GhostType & ghost_type) {
Array<UInt> & el_filter = element_filter(type, ghost_type);
el_filter.push_back(element);
return el_filter.getSize()-1;
}
/* -------------------------------------------------------------------------- */
inline UInt Material::getTangentStiffnessVoigtSize(UInt dim) const {
return (dim * (dim - 1) / 2 + dim);
}
/* -------------------------------------------------------------------------- */
inline UInt Material::getCauchyStressMatrixSize(UInt dim) const {
return (dim * dim);
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
inline void Material::gradUToF(const Matrix<Real> & grad_u,
- Matrix<Real> & F) const {
+ Matrix<Real> & F) {
AKANTU_DEBUG_ASSERT(F.size() >= grad_u.size() && grad_u.size() == dim*dim,
"The dimension of the tensor F should be greater or equal to the dimension of the tensor grad_u.");
F.eye();
for (UInt i = 0; i < dim; ++i)
for (UInt j = 0; j < dim; ++j)
F(i, j) += grad_u(i, j);
}
/* -------------------------------------------------------------------------- */
template<UInt dim >
inline void Material::computeCauchyStressOnQuad(const Matrix<Real> & F,
const Matrix<Real> & piola,
Matrix<Real> & sigma,
const Real & C33 ) const {
Real J = F.det() * sqrt(C33);
Matrix<Real> F_S(dim, dim);
F_S.mul<false, false>(F, piola);
Real constant = J ? 1./J : 0;
sigma.mul<false, true>(F_S, F, constant);
}
/* -------------------------------------------------------------------------- */
inline void Material::rightCauchy(const Matrix<Real> & F,
- Matrix<Real> & C) const {
+ Matrix<Real> & C) {
C.mul<true, false>(F, F);
}
/* -------------------------------------------------------------------------- */
inline void Material::leftCauchy(const Matrix<Real> & F,
- Matrix<Real> & B) const {
+ Matrix<Real> & B) {
B.mul<false, true>(F, F);
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
inline void Material::gradUToEpsilon(const Matrix<Real> & grad_u,
- Matrix<Real> & epsilon) const {
+ Matrix<Real> & epsilon) {
for (UInt i = 0; i < dim; ++i)
for (UInt j = 0; j < dim; ++j)
epsilon(i, j) = 0.5*(grad_u(i, j) + grad_u(j, i));
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
inline void Material::gradUToGreenStrain(const Matrix<Real> & grad_u,
- Matrix<Real> & epsilon) const {
+ Matrix<Real> & epsilon) {
epsilon.mul<true, false>(grad_u, grad_u, .5);
for (UInt i = 0; i < dim; ++i)
for (UInt j = 0; j < dim; ++j)
epsilon(i, j) += 0.5 * (grad_u(i, j) + grad_u(j, i));
}
+/* -------------------------------------------------------------------------- */
+inline Real Material::stressToVonMises(const Matrix<Real> & stress) {
+ // compute deviatoric stress
+ UInt dim = stress.cols();
+ Matrix<Real> deviatoric_stress = Matrix<Real>::eye(dim, -1. * stress.trace() / 3.);
+
+ for (UInt i = 0; i < dim; ++i)
+ for (UInt j = 0; j < dim; ++j)
+ deviatoric_stress(i,j) += stress(i,j);
+
+ // return Von Mises stress
+ return std::sqrt(3. * deviatoric_stress.doubleDot(deviatoric_stress) / 2.);
+}
+
/* ---------------------------------------------------------------------------*/
template<UInt dim>
inline void Material::SetCauchyStressArray(const Matrix<Real> & S_t, Matrix<Real> & Stress_vect) {
AKANTU_DEBUG_IN();
Stress_vect.clear();
//UInt cauchy_matrix_size = getCauchyStressArraySize(dim);
//see Finite ekement formulations for large deformation dynamic analysis, Bathe et al. IJNME vol 9, 1975, page 364 ^t\tau
/*
* 1d: [ s11 ]'
* 2d: [ s11 s22 s12 ]'
* 3d: [ s11 s22 s33 s23 s13 s12 ]
*/
for (UInt i = 0; i < dim; ++i)//diagonal terms
Stress_vect(i, 0) = S_t(i, i);
for (UInt i = 1; i < dim; ++i)// term s12 in 2D and terms s23 s13 in 3D
Stress_vect(dim+i-1, 0) = S_t(dim-i-1, dim-1);
for (UInt i = 2; i < dim; ++i)//term s13 in 3D
Stress_vect(dim+i, 0) = S_t(0, 1);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
inline void Material::setCauchyStressMatrix(const Matrix<Real> & S_t, Matrix<Real> & Stress_matrix) {
AKANTU_DEBUG_IN();
Stress_matrix.clear();
/// see Finite ekement formulations for large deformation dynamic analysis,
/// Bathe et al. IJNME vol 9, 1975, page 364 ^t\tau
for (UInt i = 0; i < dim; ++i) {
for (UInt m = 0; m < dim; ++m) {
for (UInt n = 0; n < dim; ++n) {
Stress_matrix(i * dim + m, i * dim + n) = S_t(m, n);
}
}
}
//other terms from the diagonal
/*for (UInt i = 0; i < 3 - dim; ++i) {
Stress_matrix(dim * dim + i, dim * dim + i) = S_t(dim + i, dim + i);
}*/
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
inline Element Material::convertToLocalElement(const Element & global_element) const {
UInt ge = global_element.element;
#ifndef AKANTU_NDEBUG
UInt model_mat_index = this->model->getMaterialByElement(global_element.type,
global_element.ghost_type)(ge);
UInt mat_index = this->model->getMaterialIndex(this->name);
AKANTU_DEBUG_ASSERT(model_mat_index == mat_index,
"Conversion of a global element in a local element for the wrong material "
<< this->name << std::endl);
#endif
UInt le = this->model->getMaterialLocalNumbering(global_element.type,
global_element.ghost_type)(ge);
Element tmp_quad(global_element.type,
le,
global_element.ghost_type);
return tmp_quad;
}
/* -------------------------------------------------------------------------- */
inline Element Material::convertToGlobalElement(const Element & local_element) const {
UInt le = local_element.element;
UInt ge = this->element_filter(local_element.type, local_element.ghost_type)(le);
Element tmp_quad(local_element.type,
ge,
local_element.ghost_type);
return tmp_quad;
}
/* -------------------------------------------------------------------------- */
inline QuadraturePoint Material::convertToLocalPoint(const QuadraturePoint & global_point) const {
const FEEngine & fem = this->model->getFEEngine();
UInt nb_quad = fem.getNbQuadraturePoints(global_point.type);
Element el = this->convertToLocalElement(static_cast<const Element &>(global_point));
QuadraturePoint tmp_quad(el, global_point.num_point, nb_quad);
return tmp_quad;
}
/* -------------------------------------------------------------------------- */
inline QuadraturePoint Material::convertToGlobalPoint(const QuadraturePoint & local_point) const {
const FEEngine & fem = this->model->getFEEngine();
UInt nb_quad = fem.getNbQuadraturePoints(local_point.type);
Element el = this->convertToGlobalElement(static_cast<const Element &>(local_point));
QuadraturePoint tmp_quad(el, local_point.num_point, nb_quad);
return tmp_quad;
}
/* -------------------------------------------------------------------------- */
template<ElementType type>
inline void Material::buildElementalFieldInterpolationCoodinates(__attribute__((unused)) const Matrix<Real> & coordinates,
__attribute__((unused)) Matrix<Real> & coordMatrix) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/* -------------------------------------------------------------------------- */
inline void Material::buildElementalFieldInterpolationCoodinatesLinear(const Matrix<Real> & coordinates,
Matrix<Real> & coordMatrix) {
for (UInt i = 0; i < coordinates.cols(); ++i)
coordMatrix(i, 0) = 1;
}
/* -------------------------------------------------------------------------- */
inline void Material::buildElementalFieldInterpolationCoodinatesQuadratic(const Matrix<Real> & coordinates,
Matrix<Real> & coordMatrix) {
UInt nb_quadrature_points = coordMatrix.cols();
for (UInt i = 0; i < coordinates.cols(); ++i) {
coordMatrix(i, 0) = 1;
for (UInt j = 1; j < nb_quadrature_points; ++j)
coordMatrix(i, j) = coordinates(j-1, i);
}
}
/* -------------------------------------------------------------------------- */
template<>
inline void Material::buildElementalFieldInterpolationCoodinates<_segment_2>(const Matrix<Real> & coordinates,
Matrix<Real> & coordMatrix) {
buildElementalFieldInterpolationCoodinatesLinear(coordinates, coordMatrix);
}
/* -------------------------------------------------------------------------- */
template<>
inline void Material::buildElementalFieldInterpolationCoodinates<_segment_3>(const Matrix<Real> & coordinates,
Matrix<Real> & coordMatrix) {
buildElementalFieldInterpolationCoodinatesQuadratic(coordinates, coordMatrix);
}
/* -------------------------------------------------------------------------- */
template<>
inline void Material::buildElementalFieldInterpolationCoodinates<_triangle_3>(const Matrix<Real> & coordinates,
Matrix<Real> & coordMatrix) {
buildElementalFieldInterpolationCoodinatesLinear(coordinates, coordMatrix);
}
/* -------------------------------------------------------------------------- */
template<>
inline void Material::buildElementalFieldInterpolationCoodinates<_triangle_6>(const Matrix<Real> & coordinates,
Matrix<Real> & coordMatrix) {
buildElementalFieldInterpolationCoodinatesQuadratic(coordinates, coordMatrix);
}
/* -------------------------------------------------------------------------- */
template<>
inline void Material::buildElementalFieldInterpolationCoodinates<_tetrahedron_4>(const Matrix<Real> & coordinates,
Matrix<Real> & coordMatrix) {
buildElementalFieldInterpolationCoodinatesLinear(coordinates, coordMatrix);
}
/* -------------------------------------------------------------------------- */
template<>
inline void Material::buildElementalFieldInterpolationCoodinates<_tetrahedron_10>(const Matrix<Real> & coordinates,
Matrix<Real> & coordMatrix) {
buildElementalFieldInterpolationCoodinatesQuadratic(coordinates, coordMatrix);
}
/**
* @todo Write a more efficient interpolation for quadrangles by
* dropping unnecessary quadrature points
*
*/
/* -------------------------------------------------------------------------- */
template<>
inline void Material::buildElementalFieldInterpolationCoodinates<_quadrangle_4>(const Matrix<Real> & coordinates,
Matrix<Real> & coordMatrix) {
for (UInt i = 0; i < coordinates.cols(); ++i) {
Real x = coordinates(0, i);
Real y = coordinates(1, i);
coordMatrix(i, 0) = 1;
coordMatrix(i, 1) = x;
coordMatrix(i, 2) = y;
coordMatrix(i, 3) = x * y;
}
}
/* -------------------------------------------------------------------------- */
template<>
inline void Material::buildElementalFieldInterpolationCoodinates<_quadrangle_8>(const Matrix<Real> & coordinates,
Matrix<Real> & coordMatrix) {
for (UInt i = 0; i < coordinates.cols(); ++i) {
UInt j = 0;
Real x = coordinates(0, i);
Real y = coordinates(1, i);
for (UInt e = 0; e <= 2; ++e) {
for (UInt n = 0; n <= 2; ++n) {
coordMatrix(i, j) = std::pow(x, e) * std::pow(y, n);
++j;
}
}
}
}
/* -------------------------------------------------------------------------- */
inline UInt Material::getNbDataForElements(const Array<Element> & elements,
SynchronizationTag tag) const {
if(tag == _gst_smm_stress) {
return (this->isFiniteDeformation() ? 3 : 1) * spatial_dimension * spatial_dimension *
sizeof(Real) * this->getModel().getNbQuadraturePoints(elements);
}
return 0;
}
/* -------------------------------------------------------------------------- */
inline void Material::packElementData(CommunicationBuffer & buffer,
const Array<Element> & elements,
SynchronizationTag tag) const {
if(tag == _gst_smm_stress) {
if(this->isFiniteDeformation()) {
packElementDataHelper(piola_kirchhoff_2, buffer, elements);
packElementDataHelper(gradu, buffer, elements);
}
packElementDataHelper(stress, buffer, elements);
}
}
/* -------------------------------------------------------------------------- */
inline void Material::unpackElementData(CommunicationBuffer & buffer,
const Array<Element> & elements,
SynchronizationTag tag) {
if(tag == _gst_smm_stress) {
if(this->isFiniteDeformation()) {
unpackElementDataHelper(piola_kirchhoff_2, buffer, elements);
unpackElementDataHelper(gradu, buffer, elements);
}
unpackElementDataHelper(stress, buffer, elements);
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline const T & Material::getParam(const ID & param) const {
try {
return get<T>(param);
} catch (...) {
AKANTU_EXCEPTION("No parameter " << param << " in the material " << getID());
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline void Material::setParam(const ID & param, T value) {
try {
set<T>(param, value);
} catch(...) {
AKANTU_EXCEPTION("No parameter " << param << " in the material " << getID());
}
updateInternalParameters();
}
/* -------------------------------------------------------------------------- */
template<typename T>
inline void Material::packElementDataHelper(const ElementTypeMapArray<T> & data_to_pack,
CommunicationBuffer & buffer,
const Array<Element> & elements,
const ID & fem_id) const {
DataAccessor::packElementalDataHelper<T>(data_to_pack, buffer, elements, true,
model->getFEEngine(fem_id));
}
/* -------------------------------------------------------------------------- */
template<typename T>
inline void Material::unpackElementDataHelper(ElementTypeMapArray<T> & data_to_unpack,
CommunicationBuffer & buffer,
const Array<Element> & elements,
const ID & fem_id) {
DataAccessor::unpackElementalDataHelper<T>(data_to_unpack, buffer, elements, true,
model->getFEEngine(fem_id));
}
/* -------------------------------------------------------------------------- */
template<> inline void Material::registerInternal<Real>(InternalField<Real> & vect) {
internal_vectors_real[vect.getID()] = &vect;
}
template<> inline void Material::registerInternal<UInt>(InternalField<UInt> & vect) {
internal_vectors_uint[vect.getID()] = &vect;
}
/* -------------------------------------------------------------------------- */
template<> inline void Material::unregisterInternal<Real>(InternalField<Real> & vect) {
internal_vectors_real.erase(vect.getID());
}
template<> inline void Material::unregisterInternal<UInt>(InternalField<UInt> & vect) {
internal_vectors_uint.erase(vect.getID());
}
/* -------------------------------------------------------------------------- */
inline bool Material::isInternal(const ID & id, const ElementKind & element_kind) const {
std::map<ID, InternalField<Real> *>::const_iterator internal_array =
internal_vectors_real.find(this->getID()+":"+id);
if (internal_array == internal_vectors_real.end()) return false;
if (internal_array->second->getElementKind() != element_kind) return false;
return true;
}
/* -------------------------------------------------------------------------- */
inline ElementTypeMap<UInt> Material::getInternalDataPerElem(const ID & id, const ElementKind & element_kind) const {
std::map<ID, InternalField<Real> *>::const_iterator internal_array =
internal_vectors_real.find(this->getID()+":"+id);
if (internal_array == internal_vectors_real.end()) AKANTU_EXCEPTION("cannot find internal " << id);
if (internal_array->second->getElementKind() != element_kind) AKANTU_EXCEPTION("cannot find internal " << id);
InternalField<Real> & internal = *internal_array->second;
InternalField<Real>::type_iterator it = internal.firstType(spatial_dimension, _not_ghost,element_kind);
InternalField<Real>::type_iterator last_type = internal.lastType(spatial_dimension, _not_ghost,element_kind);
ElementTypeMap<UInt> res;
for(; it != last_type; ++it) {
UInt nb_quadrature_points = 0;
if (element_kind == _ek_regular)
nb_quadrature_points = model->getFEEngine().getNbQuadraturePoints(*it);
#if defined(AKANTU_COHESIVE_ELEMENT)
else if (element_kind == _ek_cohesive)
nb_quadrature_points = model->getFEEngine("CohesiveFEEngine").getNbQuadraturePoints(*it);
#endif
res(*it) = internal.getNbComponent() * nb_quadrature_points;
}
return res;
}
diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc
index c15522ea8..f93dc12d7 100644
--- a/src/model/solid_mechanics/solid_mechanics_model.cc
+++ b/src/model/solid_mechanics/solid_mechanics_model.cc
@@ -1,1801 +1,1815 @@
/**
* @file solid_mechanics_model.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author David Simon Kammer <david.kammer@epfl.ch>
* @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Tue Jul 27 2010
* @date last modification: Fri Sep 19 2014
*
* @brief Implementation of the SolidMechanicsModel class
*
* @section LICENSE
*
* Copyright (©) 2010-2012, 2014 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_math.hh"
#include "aka_common.hh"
#include "solid_mechanics_model.hh"
#include "group_manager_inline_impl.cc"
#include "dumpable_inline_impl.hh"
#include "integration_scheme_2nd_order.hh"
#include "element_group.hh"
#include "static_communicator.hh"
#include "dof_synchronizer.hh"
#include "element_group.hh"
#include <cmath>
#ifdef AKANTU_USE_MUMPS
#include "solver_mumps.hh"
#endif
#ifdef AKANTU_USE_PETSC
#include "solver_petsc.hh"
#include "petsc_matrix.hh"
#endif
#ifdef AKANTU_USE_IOHELPER
# include "dumper_field.hh"
# include "dumper_paraview.hh"
# include "dumper_homogenizing_field.hh"
# include "dumper_material_internal_field.hh"
# include "dumper_elemental_field.hh"
# include "dumper_material_padders.hh"
# include "dumper_element_partition.hh"
# include "dumper_iohelper.hh"
#endif
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
const SolidMechanicsModelOptions default_solid_mechanics_model_options(_explicit_lumped_mass, false);
/* -------------------------------------------------------------------------- */
/**
* A solid mechanics model need a mesh and a dimension to be created. the model
* by it self can not do a lot, the good init functions should be called in
* order to configure the model depending on what we want to do.
*
* @param mesh mesh representing the model we want to simulate
* @param dim spatial dimension of the problem, if dim = 0 (default value) the
* dimension of the problem is assumed to be the on of the mesh
* @param id an id to identify the model
*/
SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh,
UInt dim,
const ID & id,
const MemoryID & memory_id) :
Model(mesh, dim, id, memory_id),
BoundaryCondition<SolidMechanicsModel>(),
time_step(NAN), f_m2a(1.0),
mass_matrix(NULL),
velocity_damping_matrix(NULL),
stiffness_matrix(NULL),
jacobian_matrix(NULL),
material_index("material index", id),
material_local_numbering("material local numbering", id),
material_selector(new DefaultMaterialSelector(material_index)),
is_default_material_selector(true),
integrator(NULL),
increment_flag(false), solver(NULL),
synch_parallel(NULL),
are_materials_instantiated(false) {
AKANTU_DEBUG_IN();
createSynchronizerRegistry(this);
registerFEEngineObject<MyFEEngineType>("SolidMechanicsFEEngine", mesh, spatial_dimension);
this->displacement = NULL;
this->mass = NULL;
this->velocity = NULL;
this->acceleration = NULL;
this->force = NULL;
this->residual = NULL;
this->blocked_dofs = NULL;
this->increment = NULL;
this->increment_acceleration = NULL;
this->dof_synchronizer = NULL;
this->previous_displacement = NULL;
materials.clear();
mesh.registerEventHandler(*this);
#if defined(AKANTU_USE_IOHELPER)
this->mesh.registerDumper<DumperParaview>("paraview_all", id, true);
this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_regular);
#endif
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
SolidMechanicsModel::~SolidMechanicsModel() {
AKANTU_DEBUG_IN();
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
delete *mat_it;
}
materials.clear();
delete integrator;
delete solver;
delete mass_matrix;
delete velocity_damping_matrix;
if(stiffness_matrix && stiffness_matrix != jacobian_matrix)
delete stiffness_matrix;
delete jacobian_matrix;
delete synch_parallel;
if(is_default_material_selector) {
delete material_selector;
material_selector = NULL;
}
AKANTU_DEBUG_OUT();
}
void SolidMechanicsModel::setTimeStep(Real time_step) {
this->time_step = time_step;
#if defined(AKANTU_USE_IOHELPER)
this->mesh.getDumper().setTimeStep(time_step);
#endif
}
/* -------------------------------------------------------------------------- */
/* Initialisation */
/* -------------------------------------------------------------------------- */
/**
* This function groups many of the initialization in on function. For most of
* basics case the function should be enough. The functions initialize the
* model, the internal vectors, set them to 0, and depending on the parameters
* it also initialize the explicit or implicit solver.
*
* @param material_file the file containing the materials to use
* @param method the analysis method wanted. See the akantu::AnalysisMethod for
* the different possibilities
*/
void SolidMechanicsModel::initFull(const ModelOptions & options) {
Model::initFull(options);
const SolidMechanicsModelOptions & smm_options =
dynamic_cast<const SolidMechanicsModelOptions &>(options);
method = smm_options.analysis_method;
// initialize the vectors
initArrays();
// set the initial condition to 0
force->clear();
velocity->clear();
acceleration->clear();
displacement->clear();
// initialize pcb
if(pbc_pair.size()!=0)
initPBC();
// initialize the time integration schemes
switch(method) {
case _explicit_lumped_mass:
initExplicit();
break;
case _explicit_consistent_mass:
initSolver();
initExplicit();
break;
case _implicit_dynamic:
initImplicit(true);
break;
case _static:
initImplicit(false);
break;
default:
AKANTU_EXCEPTION("analysis method not recognised by SolidMechanicsModel");
break;
}
// initialize the materials
if(this->parser->getLastParsedFile() != "") {
instantiateMaterials();
}
if(!smm_options.no_init_materials) {
initMaterials();
}
if(increment_flag)
initBC(*this, *displacement, *increment, *force);
else
initBC(*this, *displacement, *force);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initParallel(MeshPartition * partition,
DataAccessor * data_accessor) {
AKANTU_DEBUG_IN();
if (data_accessor == NULL) data_accessor = this;
synch_parallel = &createParallelSynch(partition,data_accessor);
synch_registry->registerSynchronizer(*synch_parallel, _gst_material_id);
synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_mass);
synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_stress);
synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_boundary);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initFEEngineBoundary() {
FEEngine & fem_boundary = getFEEngineBoundary();
fem_boundary.initShapeFunctions(_not_ghost);
fem_boundary.initShapeFunctions(_ghost);
fem_boundary.computeNormalsOnControlPoints(_not_ghost);
fem_boundary.computeNormalsOnControlPoints(_ghost);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initExplicit(AnalysisMethod analysis_method) {
AKANTU_DEBUG_IN();
//in case of switch from implicit to explicit
if(!this->isExplicit())
method = analysis_method;
if (integrator) delete integrator;
integrator = new CentralDifference();
UInt nb_nodes = acceleration->getSize();
UInt nb_degree_of_freedom = acceleration->getNbComponent();
std::stringstream sstr; sstr << id << ":increment_acceleration";
increment_acceleration = &(alloc<Real>(sstr.str(), nb_nodes, nb_degree_of_freedom, Real()));
AKANTU_DEBUG_OUT();
}
void SolidMechanicsModel::initArraysPreviousDisplacment() {
AKANTU_DEBUG_IN();
SolidMechanicsModel::setIncrementFlagOn();
UInt nb_nodes = mesh.getNbNodes();
std::stringstream sstr_disp_t;
sstr_disp_t << id << ":previous_displacement";
previous_displacement = &(alloc<Real > (sstr_disp_t.str(), nb_nodes, spatial_dimension, 0.));
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Allocate all the needed vectors. By default their are not necessarily set to
* 0
*
*/
void SolidMechanicsModel::initArrays() {
AKANTU_DEBUG_IN();
UInt nb_nodes = mesh.getNbNodes();
std::stringstream sstr_disp; sstr_disp << id << ":displacement";
// std::stringstream sstr_mass; sstr_mass << id << ":mass";
std::stringstream sstr_velo; sstr_velo << id << ":velocity";
std::stringstream sstr_acce; sstr_acce << id << ":acceleration";
std::stringstream sstr_forc; sstr_forc << id << ":force";
std::stringstream sstr_resi; sstr_resi << id << ":residual";
std::stringstream sstr_boun; sstr_boun << id << ":blocked_dofs";
displacement = &(alloc<Real>(sstr_disp.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE));
// mass = &(alloc<Real>(sstr_mass.str(), nb_nodes, spatial_dimension, 0));
velocity = &(alloc<Real>(sstr_velo.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE));
acceleration = &(alloc<Real>(sstr_acce.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE));
force = &(alloc<Real>(sstr_forc.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE));
residual = &(alloc<Real>(sstr_resi.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE));
blocked_dofs = &(alloc<bool>(sstr_boun.str(), nb_nodes, spatial_dimension, false));
std::stringstream sstr_curp; sstr_curp << id << ":current_position";
current_position = &(alloc<Real>(sstr_curp.str(), 0, spatial_dimension, REAL_INIT_VALUE));
for(UInt g = _not_ghost; g <= _ghost; ++g) {
GhostType gt = (GhostType) g;
Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_not_defined);
Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_not_defined);
for(; it != end; ++it) {
UInt nb_element = mesh.getNbElement(*it, gt);
material_index.alloc(nb_element, 1, *it, gt);
material_local_numbering.alloc(nb_element, 1, *it, gt);
}
}
dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension);
dof_synchronizer->initLocalDOFEquationNumbers();
dof_synchronizer->initGlobalDOFEquationNumbers();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Initialize the model,basically it pre-compute the shapes, shapes derivatives
* and jacobian
*
*/
void SolidMechanicsModel::initModel() {
/// \todo add the current position as a parameter to initShapeFunctions for
/// large deformation
getFEEngine().initShapeFunctions(_not_ghost);
getFEEngine().initShapeFunctions(_ghost);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initPBC() {
Model::initPBC();
registerPBCSynchronizer();
// as long as there are ones on the diagonal of the matrix, we can put boudandary true for slaves
std::map<UInt, UInt>::iterator it = pbc_pair.begin();
std::map<UInt, UInt>::iterator end = pbc_pair.end();
UInt dim = mesh.getSpatialDimension();
while(it != end) {
for (UInt i=0; i<dim; ++i)
(*blocked_dofs)((*it).first,i) = true;
++it;
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::registerPBCSynchronizer(){
PBCSynchronizer * synch = new PBCSynchronizer(pbc_pair);
synch_registry->registerSynchronizer(*synch, _gst_smm_uv);
synch_registry->registerSynchronizer(*synch, _gst_smm_mass);
synch_registry->registerSynchronizer(*synch, _gst_smm_res);
synch_registry->registerSynchronizer(*synch, _gst_for_dump);
changeLocalEquationNumberForPBC(pbc_pair, mesh.getSpatialDimension());
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateCurrentPosition() {
AKANTU_DEBUG_IN();
UInt nb_nodes = mesh.getNbNodes();
current_position->resize(nb_nodes);
Real * current_position_val = current_position->storage();
Real * position_val = mesh.getNodes().storage();
Real * displacement_val = displacement->storage();
/// compute current_position = initial_position + displacement
memcpy(current_position_val, position_val, nb_nodes*spatial_dimension*sizeof(Real));
for (UInt n = 0; n < nb_nodes*spatial_dimension; ++n) {
*current_position_val++ += *displacement_val++;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initializeUpdateResidualData() {
AKANTU_DEBUG_IN();
UInt nb_nodes = mesh.getNbNodes();
residual->resize(nb_nodes);
/// copy the forces in residual for boundary conditions
memcpy(residual->storage(), force->storage(), nb_nodes*spatial_dimension*sizeof(Real));
// start synchronization
synch_registry->asynchronousSynchronize(_gst_smm_uv);
synch_registry->waitEndSynchronize(_gst_smm_uv);
updateCurrentPosition();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/* Explicit scheme */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/**
* This function compute the second member of the motion equation. That is to
* say the sum of forces @f$ r = F_{ext} - F_{int} @f$. @f$ F_{ext} @f$ is given
* by the user in the force vector, and @f$ F_{int} @f$ is computed as @f$
* F_{int} = \int_{\Omega} N \sigma d\Omega@f$
*
*/
void SolidMechanicsModel::updateResidual(bool need_initialize) {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Assemble the internal forces");
// f = f_ext - f_int
// f = f_ext
if(need_initialize) initializeUpdateResidualData();
AKANTU_DEBUG_INFO("Compute local stresses");
std::vector<Material *>::iterator mat_it;
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.computeAllStresses(_not_ghost);
}
#ifdef AKANTU_DAMAGE_NON_LOCAL
/* ------------------------------------------------------------------------ */
/* Computation of the non local part */
synch_registry->asynchronousSynchronize(_gst_mnl_for_average);
AKANTU_DEBUG_INFO("Compute non local stresses for local elements");
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.computeAllNonLocalStresses(_not_ghost);
}
AKANTU_DEBUG_INFO("Wait distant non local stresses");
synch_registry->waitEndSynchronize(_gst_mnl_for_average);
AKANTU_DEBUG_INFO("Compute non local stresses for ghosts elements");
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.computeAllNonLocalStresses(_ghost);
}
#endif
/* ------------------------------------------------------------------------ */
/* assembling the forces internal */
// communicate the stress
AKANTU_DEBUG_INFO("Send data for residual assembly");
synch_registry->asynchronousSynchronize(_gst_smm_stress);
AKANTU_DEBUG_INFO("Assemble residual for local elements");
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.assembleResidual(_not_ghost);
}
AKANTU_DEBUG_INFO("Wait distant stresses");
// finalize communications
synch_registry->waitEndSynchronize(_gst_smm_stress);
AKANTU_DEBUG_INFO("Assemble residual for ghost elements");
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.assembleResidual(_ghost);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::computeStresses() {
if (isExplicit()) {
// start synchronization
synch_registry->asynchronousSynchronize(_gst_smm_uv);
synch_registry->waitEndSynchronize(_gst_smm_uv);
// compute stresses on all local elements for each materials
std::vector<Material *>::iterator mat_it;
for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.computeAllStresses(_not_ghost);
}
/* ------------------------------------------------------------------------ */
#ifdef AKANTU_DAMAGE_NON_LOCAL
/* Computation of the non local part */
synch_registry->asynchronousSynchronize(_gst_mnl_for_average);
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.computeAllNonLocalStresses(_not_ghost);
}
synch_registry->waitEndSynchronize(_gst_mnl_for_average);
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.computeAllNonLocalStresses(_ghost);
}
#endif
} else {
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
mat.computeAllStressesFromTangentModuli(_not_ghost);
}
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateResidualInternal() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Update the residual");
// f = f_ext - f_int - Ma - Cv = r - Ma - Cv;
if(method != _static) {
// f -= Ma
if(mass_matrix) {
// if full mass_matrix
Array<Real> * Ma = new Array<Real>(*acceleration, true, "Ma");
*Ma *= *mass_matrix;
/// \todo check unit conversion for implicit dynamics
// *Ma /= f_m2a
*residual -= *Ma;
delete Ma;
} else if (mass) {
// else lumped mass
UInt nb_nodes = acceleration->getSize();
UInt nb_degree_of_freedom = acceleration->getNbComponent();
Real * mass_val = mass->storage();
Real * accel_val = acceleration->storage();
Real * res_val = residual->storage();
bool * blocked_dofs_val = blocked_dofs->storage();
for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) {
if(!(*blocked_dofs_val)) {
*res_val -= *accel_val * *mass_val /f_m2a;
}
blocked_dofs_val++;
res_val++;
mass_val++;
accel_val++;
}
} else {
AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix.");
}
// f -= Cv
if(velocity_damping_matrix) {
Array<Real> * Cv = new Array<Real>(*velocity);
*Cv *= *velocity_damping_matrix;
*residual -= *Cv;
delete Cv;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateAcceleration() {
AKANTU_DEBUG_IN();
updateResidualInternal();
if(method == _explicit_lumped_mass) {
/* residual = residual_{n+1} - M * acceleration_n therefore
solution = increment acceleration not acceleration */
solveLumped(*increment_acceleration,
*mass,
*residual,
*blocked_dofs,
f_m2a);
} else if (method == _explicit_consistent_mass) {
solve<NewmarkBeta::_acceleration_corrector>(*increment_acceleration);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::solveLumped(Array<Real> & x,
const Array<Real> & A,
const Array<Real> & b,
const Array<bool> & blocked_dofs,
Real alpha) {
Real * A_val = A.storage();
Real * b_val = b.storage();
Real * x_val = x.storage();
bool * blocked_dofs_val = blocked_dofs.storage();
UInt nb_degrees_of_freedom = x.getSize() * x.getNbComponent();
for (UInt n = 0; n < nb_degrees_of_freedom; ++n) {
if(!(*blocked_dofs_val)) {
*x_val = alpha * (*b_val / *A_val);
}
x_val++;
A_val++;
b_val++;
blocked_dofs_val++;
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::explicitPred() {
AKANTU_DEBUG_IN();
if(increment_flag) {
if(previous_displacement) increment->copy(*previous_displacement);
else increment->copy(*displacement);
}
AKANTU_DEBUG_ASSERT(integrator,"itegrator should have been allocated: "
<< "have called initExplicit ? "
<< "or initImplicit ?");
integrator->integrationSchemePred(time_step,
*displacement,
*velocity,
*acceleration,
*blocked_dofs);
if(increment_flag) {
Real * inc_val = increment->storage();
Real * dis_val = displacement->storage();
UInt nb_degree_of_freedom = displacement->getSize() * displacement->getNbComponent();
for (UInt n = 0; n < nb_degree_of_freedom; ++n) {
*inc_val = *dis_val - *inc_val;
inc_val++;
dis_val++;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::explicitCorr() {
AKANTU_DEBUG_IN();
integrator->integrationSchemeCorrAccel(time_step,
*displacement,
*velocity,
*acceleration,
*blocked_dofs,
*increment_acceleration);
if(previous_displacement) previous_displacement->copy(*displacement);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::solveStep() {
AKANTU_DEBUG_IN();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeSolveStepEvent(method));
this->explicitPred();
this->updateResidual();
this->updateAcceleration();
this->explicitCorr();
EventManager::sendEvent(SolidMechanicsModelEvent::AfterSolveStepEvent(method));
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/* Implicit scheme */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/**
* Initialize the solver and create the sparse matrices needed.
*
*/
void SolidMechanicsModel::initSolver(__attribute__((unused)) SolverOptions & options) {
#if !defined(AKANTU_USE_MUMPS) && !defined(AKANTU_USE_PETSC)// or other solver in the future \todo add AKANTU_HAS_SOLVER in CMake
AKANTU_DEBUG_ERROR("You should at least activate one solver.");
#else
UInt nb_global_nodes = mesh.getNbGlobalNodes();
delete jacobian_matrix;
std::stringstream sstr; sstr << id << ":jacobian_matrix";
#ifdef AKANTU_USE_PETSC
jacobian_matrix = new PETScMatrix(nb_global_nodes * spatial_dimension, _symmetric, sstr.str(), memory_id);
#else
jacobian_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension, _symmetric, sstr.str(), memory_id);
#endif //AKANTU_USE PETSC
jacobian_matrix->buildProfile(mesh, *dof_synchronizer, spatial_dimension);
if (!isExplicit()) {
delete stiffness_matrix;
std::stringstream sstr_sti; sstr_sti << id << ":stiffness_matrix";
#ifdef AKANTU_USE_PETSC
stiffness_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension, _symmetric, sstr.str(), memory_id);
stiffness_matrix->buildProfile(mesh, *dof_synchronizer, spatial_dimension);
#else
stiffness_matrix = new SparseMatrix(*jacobian_matrix, sstr_sti.str(), memory_id);
#endif //AKANTU_USE_PETSC
}
delete solver;
std::stringstream sstr_solv; sstr_solv << id << ":solver";
#ifdef AKANTU_USE_PETSC
solver = new SolverPETSc(*jacobian_matrix, sstr_solv.str());
#elif defined(AKANTU_USE_MUMPS)
solver = new SolverMumps(*jacobian_matrix, sstr_solv.str());
dof_synchronizer->initScatterGatherCommunicationScheme();
#else
AKANTU_DEBUG_ERROR("You should at least activate one solver.");
#endif //AKANTU_USE_MUMPS
if(solver)
solver->initialize(options);
#endif //AKANTU_HAS_SOLVER
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initJacobianMatrix() {
#if defined(AKANTU_USE_MUMPS) && !defined(AKANTU_USE_PETSC)
// @todo make it more flexible: this is an ugly patch to treat the case of non
// fix profile of the K matrix
delete jacobian_matrix;
std::stringstream sstr_sti; sstr_sti << id << ":jacobian_matrix";
jacobian_matrix = new SparseMatrix(*stiffness_matrix, sstr_sti.str(), memory_id);
std::stringstream sstr_solv; sstr_solv << id << ":solver";
delete solver;
solver = new SolverMumps(*jacobian_matrix, sstr_solv.str());
if(solver)
solver->initialize(_solver_no_options);
#else
AKANTU_DEBUG_ERROR("You need to activate the solver mumps.");
#endif
}
/* -------------------------------------------------------------------------- */
/**
* Initialize the implicit solver, either for dynamic or static cases,
*
* @param dynamic
*/
void SolidMechanicsModel::initImplicit(bool dynamic, SolverOptions & solver_options) {
AKANTU_DEBUG_IN();
method = dynamic ? _implicit_dynamic : _static;
if (!increment) setIncrementFlagOn();
initSolver(solver_options);
if(method == _implicit_dynamic) {
if(integrator) delete integrator;
integrator = new TrapezoidalRule2();
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::initialAcceleration() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Solving Ma = f");
Solver * acc_solver = NULL;
std::stringstream sstr; sstr << id << ":tmp_mass_matrix";
SparseMatrix * tmp_mass = new SparseMatrix(*mass_matrix, sstr.str(), memory_id);
#ifdef AKANTU_USE_MUMPS
std::stringstream sstr_solver; sstr << id << ":solver_mass_matrix";
acc_solver = new SolverMumps(*mass_matrix, sstr_solver.str());
dof_synchronizer->initScatterGatherCommunicationScheme();
#else
AKANTU_DEBUG_ERROR("You should at least activate one solver.");
#endif //AKANTU_USE_MUMPS
acc_solver->initialize();
tmp_mass->applyBoundary(*blocked_dofs);
acc_solver->setRHS(*residual);
acc_solver->solve(*acceleration);
delete acc_solver;
delete tmp_mass;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::assembleStiffnessMatrix() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Assemble the new stiffness matrix.");
stiffness_matrix->clear();
// call compute stiffness matrix on each local elements
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
(*mat_it)->assembleStiffnessMatrix(_not_ghost);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
SparseMatrix & SolidMechanicsModel::initVelocityDampingMatrix() {
if(!velocity_damping_matrix)
velocity_damping_matrix =
new SparseMatrix(*jacobian_matrix, id + ":velocity_damping_matrix", memory_id);
return *velocity_damping_matrix;
}
/* -------------------------------------------------------------------------- */
template<>
bool SolidMechanicsModel::testConvergence<_scc_increment>(Real tolerance, Real & error){
AKANTU_DEBUG_IN();
UInt nb_nodes = displacement->getSize();
UInt nb_degree_of_freedom = displacement->getNbComponent();
error = 0;
Real norm[2] = {0., 0.};
Real * increment_val = increment->storage();
bool * blocked_dofs_val = blocked_dofs->storage();
Real * displacement_val = displacement->storage();
for (UInt n = 0; n < nb_nodes; ++n) {
bool is_local_node = mesh.isLocalOrMasterNode(n);
for (UInt d = 0; d < nb_degree_of_freedom; ++d) {
if(!(*blocked_dofs_val) && is_local_node) {
norm[0] += *increment_val * *increment_val;
norm[1] += *displacement_val * *displacement_val;
}
blocked_dofs_val++;
increment_val++;
displacement_val++;
}
}
StaticCommunicator::getStaticCommunicator().allReduce(norm, 2, _so_sum);
norm[0] = sqrt(norm[0]);
norm[1] = sqrt(norm[1]);
AKANTU_DEBUG_ASSERT(!Math::isnan(norm[0]), "Something goes wrong in the solve phase");
if (norm[1] < Math::getTolerance()) {
error = norm[0];
AKANTU_DEBUG_OUT();
// cout<<"Error 1: "<<error<<endl;
return error < tolerance;
}
AKANTU_DEBUG_OUT();
if(norm[1] > Math::getTolerance())
error = norm[0] / norm[1];
else
error = norm[0]; //In case the total displacement is zero!
// cout<<"Error 2: "<<error<<endl;
return (error < tolerance);
}
/* -------------------------------------------------------------------------- */
template<>
bool SolidMechanicsModel::testConvergence<_scc_residual>(Real tolerance, Real & norm) {
AKANTU_DEBUG_IN();
UInt nb_nodes = residual->getSize();
norm = 0;
Real * residual_val = residual->storage();
bool * blocked_dofs_val = blocked_dofs->storage();
for (UInt n = 0; n < nb_nodes; ++n) {
bool is_local_node = mesh.isLocalOrMasterNode(n);
if(is_local_node) {
for (UInt d = 0; d < spatial_dimension; ++d) {
if(!(*blocked_dofs_val)) {
norm += *residual_val * *residual_val;
}
blocked_dofs_val++;
residual_val++;
}
} else {
blocked_dofs_val += spatial_dimension;
residual_val += spatial_dimension;
}
}
StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum);
norm = sqrt(norm);
AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something goes wrong in the solve phase");
AKANTU_DEBUG_OUT();
return (norm < tolerance);
}
/* -------------------------------------------------------------------------- */
template<>
bool SolidMechanicsModel::testConvergence<_scc_residual_mass_wgh>(Real tolerance,
Real & norm) {
AKANTU_DEBUG_IN();
UInt nb_nodes = residual->getSize();
norm = 0;
Real * residual_val = residual->storage();
Real * mass_val = this->mass->storage();
bool * blocked_dofs_val = blocked_dofs->storage();
for (UInt n = 0; n < nb_nodes; ++n) {
bool is_local_node = mesh.isLocalOrMasterNode(n);
if(is_local_node) {
for (UInt d = 0; d < spatial_dimension; ++d) {
if(!(*blocked_dofs_val)) {
norm += *residual_val * *residual_val/(*mass_val * *mass_val);
}
blocked_dofs_val++;
residual_val++;
mass_val++;
}
} else {
blocked_dofs_val += spatial_dimension;
residual_val += spatial_dimension;
mass_val += spatial_dimension;
}
}
StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum);
norm = sqrt(norm);
AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something goes wrong in the solve phase");
AKANTU_DEBUG_OUT();
return (norm < tolerance);
}
/* -------------------------------------------------------------------------- */
bool SolidMechanicsModel::testConvergenceResidual(Real tolerance){
AKANTU_DEBUG_IN();
Real error=0;
bool res = this->testConvergence<_scc_residual>(tolerance, error);
AKANTU_DEBUG_OUT();
return res;
}
/* -------------------------------------------------------------------------- */
bool SolidMechanicsModel::testConvergenceResidual(Real tolerance, Real & error){
AKANTU_DEBUG_IN();
bool res = this->testConvergence<_scc_residual>(tolerance, error);
AKANTU_DEBUG_OUT();
return res;
}
/* -------------------------------------------------------------------------- */
bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance){
AKANTU_DEBUG_IN();
Real error=0;
bool res = this->testConvergence<_scc_increment>(tolerance, error);
AKANTU_DEBUG_OUT();
return res;
}
/* -------------------------------------------------------------------------- */
bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance, Real & error){
AKANTU_DEBUG_IN();
bool res = this->testConvergence<_scc_increment>(tolerance, error);
AKANTU_DEBUG_OUT();
return res;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::implicitPred() {
AKANTU_DEBUG_IN();
if(previous_displacement) previous_displacement->copy(*displacement);
if(method == _implicit_dynamic)
integrator->integrationSchemePred(time_step,
*displacement,
*velocity,
*acceleration,
*blocked_dofs);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::implicitCorr() {
AKANTU_DEBUG_IN();
if(method == _implicit_dynamic) {
integrator->integrationSchemeCorrDispl(time_step,
*displacement,
*velocity,
*acceleration,
*blocked_dofs,
*increment);
} else {
UInt nb_nodes = displacement->getSize();
UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes;
Real * incr_val = increment->storage();
Real * disp_val = displacement->storage();
bool * boun_val = blocked_dofs->storage();
for (UInt j = 0; j < nb_degree_of_freedom; ++j, ++disp_val, ++incr_val, ++boun_val){
*incr_val *= (1. - *boun_val);
*disp_val += *incr_val;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateIncrement() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(previous_displacement,"The previous displacement has to be initialized."
<< " Are you working with Finite or Ineslactic deformations?");
UInt nb_nodes = displacement->getSize();
UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes;
Real * incr_val = increment->storage();
Real * disp_val = displacement->storage();
Real * prev_disp_val = previous_displacement->storage();
for (UInt j = 0; j < nb_degree_of_freedom; ++j, ++disp_val, ++incr_val, ++prev_disp_val)
*incr_val = (*disp_val - *prev_disp_val);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updatePreviousDisplacement() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(previous_displacement,"The previous displacement has to be initialized."
<< " Are you working with Finite or Ineslactic deformations?");
previous_displacement->copy(*displacement);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/* Information */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::synchronizeBoundaries() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(synch_registry,"Synchronizer registry was not initialized."
<< " Did you call initParallel?");
synch_registry->synchronize(_gst_smm_boundary);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::synchronizeResidual() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(synch_registry,"Synchronizer registry was not initialized."
<< " Did you call initPBC?");
synch_registry->synchronize(_gst_smm_res);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::setIncrementFlagOn() {
AKANTU_DEBUG_IN();
if(!increment) {
UInt nb_nodes = mesh.getNbNodes();
std::stringstream sstr_inc; sstr_inc << id << ":increment";
increment = &(alloc<Real>(sstr_inc.str(), nb_nodes, spatial_dimension, 0.));
}
increment_flag = true;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getStableTimeStep() {
AKANTU_DEBUG_IN();
Real min_dt = getStableTimeStep(_not_ghost);
/// reduction min over all processors
StaticCommunicator::getStaticCommunicator().allReduce(&min_dt, 1, _so_min);
AKANTU_DEBUG_OUT();
return min_dt;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getStableTimeStep(const GhostType & ghost_type) {
AKANTU_DEBUG_IN();
Material ** mat_val = &(materials.at(0));
Real min_dt = std::numeric_limits<Real>::max();
updateCurrentPosition();
Element elem;
elem.ghost_type = ghost_type;
elem.kind = _ek_regular;
Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type);
Mesh::type_iterator end = mesh.lastType(spatial_dimension, ghost_type);
for(; it != end; ++it) {
elem.type = *it;
UInt nb_nodes_per_element = mesh.getNbNodesPerElement(*it);
UInt nb_element = mesh.getNbElement(*it);
Array<UInt>::const_scalar_iterator mat_indexes = material_index(*it, ghost_type).begin();
Array<UInt>::const_scalar_iterator mat_loc_num = material_local_numbering(*it, ghost_type).begin();
Array<Real> X(0, nb_nodes_per_element*spatial_dimension);
FEEngine::extractNodalToElementField(mesh, *current_position,
X, *it, _not_ghost);
Array<Real>::matrix_iterator X_el =
X.begin(spatial_dimension, nb_nodes_per_element);
for (UInt el = 0; el < nb_element; ++el, ++X_el, ++mat_indexes, ++mat_loc_num) {
elem.element = *mat_loc_num;
Real el_h = getFEEngine().getElementInradius(*X_el, *it);
Real el_c = mat_val[*mat_indexes]->getCelerity(elem);
Real el_dt = el_h / el_c;
min_dt = std::min(min_dt, el_dt);
}
}
AKANTU_DEBUG_OUT();
return min_dt;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getPotentialEnergy() {
AKANTU_DEBUG_IN();
Real energy = 0.;
/// call update residual on each local elements
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
energy += (*mat_it)->getPotentialEnergy();
}
/// reduction sum over all processors
StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_sum);
AKANTU_DEBUG_OUT();
return energy;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getKineticEnergy() {
AKANTU_DEBUG_IN();
if (!mass && !mass_matrix)
AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix.");
Real ekin = 0.;
UInt nb_nodes = mesh.getNbNodes();
Real * vel_val = velocity->storage();
Real * mass_val = mass->storage();
for (UInt n = 0; n < nb_nodes; ++n) {
Real mv2 = 0;
bool is_local_node = mesh.isLocalOrMasterNode(n);
bool is_not_pbc_slave_node = !isPBCSlaveNode(n);
bool count_node = is_local_node && is_not_pbc_slave_node;
for (UInt i = 0; i < spatial_dimension; ++i) {
if (count_node)
mv2 += *vel_val * *vel_val * *mass_val;
vel_val++;
mass_val++;
}
ekin += mv2;
}
StaticCommunicator::getStaticCommunicator().allReduce(&ekin, 1, _so_sum);
AKANTU_DEBUG_OUT();
return ekin * .5;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getKineticEnergy(const ElementType & type, UInt index) {
AKANTU_DEBUG_IN();
UInt nb_quadrature_points = getFEEngine().getNbQuadraturePoints(type);
Array<Real> vel_on_quad(nb_quadrature_points, spatial_dimension);
Array<UInt> filter_element(1, 1, index);
getFEEngine().interpolateOnQuadraturePoints(*velocity, vel_on_quad,
spatial_dimension,
type, _not_ghost,
filter_element);
Array<Real>::vector_iterator vit = vel_on_quad.begin(spatial_dimension);
Array<Real>::vector_iterator vend = vel_on_quad.end(spatial_dimension);
Vector<Real> rho_v2(nb_quadrature_points);
Real rho = materials[material_index(type)(index)]->getRho();
for (UInt q = 0; vit != vend; ++vit, ++q) {
rho_v2(q) = rho * vit->dot(*vit);
}
AKANTU_DEBUG_OUT();
return .5*getFEEngine().integrate(rho_v2, type, index);
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getExternalWork() {
AKANTU_DEBUG_IN();
Real * velo = velocity->storage();
Real * forc = force->storage();
Real * resi = residual->storage();
bool * boun = blocked_dofs->storage();
Real work = 0.;
UInt nb_nodes = mesh.getNbNodes();
for (UInt n = 0; n < nb_nodes; ++n) {
bool is_local_node = mesh.isLocalOrMasterNode(n);
bool is_not_pbc_slave_node = !isPBCSlaveNode(n);
bool count_node = is_local_node && is_not_pbc_slave_node;
for (UInt i = 0; i < spatial_dimension; ++i) {
if (count_node) {
if(*boun)
work -= *resi * *velo * time_step;
else
work += *forc * *velo * time_step;
}
++velo;
++forc;
++resi;
++boun;
}
}
StaticCommunicator::getStaticCommunicator().allReduce(&work, 1, _so_sum);
AKANTU_DEBUG_OUT();
return work;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getEnergy(const std::string & energy_id) {
AKANTU_DEBUG_IN();
if (energy_id == "kinetic") {
return getKineticEnergy();
} else if (energy_id == "external work"){
return getExternalWork();
}
Real energy = 0.;
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
energy += (*mat_it)->getEnergy(energy_id);
}
/// reduction sum over all processors
StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_sum);
AKANTU_DEBUG_OUT();
return energy;
}
/* -------------------------------------------------------------------------- */
Real SolidMechanicsModel::getEnergy(const std::string & energy_id,
const ElementType & type,
UInt index){
AKANTU_DEBUG_IN();
if (energy_id == "kinetic") {
return getKineticEnergy(type, index);
}
std::vector<Material *>::iterator mat_it;
UInt mat_index = this->material_index(type, _not_ghost)(index);
UInt mat_loc_num = this->material_local_numbering(type, _not_ghost)(index);
Real energy = this->materials[mat_index]->getEnergy(energy_id, type, mat_loc_num);
AKANTU_DEBUG_OUT();
return energy;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onElementsAdded(const Array<Element> & element_list,
const NewElementsEvent & event) {
AKANTU_DEBUG_IN();
this->getFEEngine().initShapeFunctions(_not_ghost);
this->getFEEngine().initShapeFunctions(_ghost);
for(UInt g = _not_ghost; g <= _ghost; ++g) {
GhostType gt = (GhostType) g;
Mesh::type_iterator it = this->mesh.firstType(spatial_dimension, gt, _ek_not_defined);
Mesh::type_iterator end = this->mesh.lastType(spatial_dimension, gt, _ek_not_defined);
for(; it != end; ++it) {
UInt nb_element = this->mesh.getNbElement(*it, gt);
if(!material_index.exists(*it, gt)) {
this->material_index .alloc(nb_element, 1, *it, gt);
this->material_local_numbering.alloc(nb_element, 1, *it, gt);
} else {
this->material_index (*it, gt).resize(nb_element);
this->material_local_numbering(*it, gt).resize(nb_element);
}
}
}
Array<Element>::const_iterator<Element> it = element_list.begin();
Array<Element>::const_iterator<Element> end = element_list.end();
ElementTypeMapArray<UInt> filter("new_element_filter", this->getID());
for (UInt el = 0; it != end; ++it, ++el) {
const Element & elem = *it;
if(!filter.exists(elem.type, elem.ghost_type))
filter.alloc(0, 1, elem.type, elem.ghost_type);
filter(elem.type, elem.ghost_type).push_back(elem.element);
}
this->assignMaterialToElements(&filter);
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
(*mat_it)->onElementsAdded(element_list, event);
}
if(method == _explicit_lumped_mass) this->assembleMassLumped();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onElementsRemoved(__attribute__((unused)) const Array<Element> & element_list,
const ElementTypeMapArray<UInt> & new_numbering,
const RemovedElementsEvent & event) {
this->getFEEngine().initShapeFunctions(_not_ghost);
this->getFEEngine().initShapeFunctions(_ghost);
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
(*mat_it)->onElementsRemoved(element_list, new_numbering, event);
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onNodesAdded(const Array<UInt> & nodes_list,
__attribute__((unused)) const NewNodesEvent & event) {
AKANTU_DEBUG_IN();
UInt nb_nodes = mesh.getNbNodes();
if(displacement) displacement->resize(nb_nodes);
if(mass ) mass ->resize(nb_nodes);
if(velocity ) velocity ->resize(nb_nodes);
if(acceleration) acceleration->resize(nb_nodes);
if(force ) force ->resize(nb_nodes);
if(residual ) residual ->resize(nb_nodes);
if(blocked_dofs) blocked_dofs->resize(nb_nodes);
if(previous_displacement) previous_displacement->resize(nb_nodes);
if(increment_acceleration) increment_acceleration->resize(nb_nodes);
if(increment) increment->resize(nb_nodes);
if(current_position) current_position->resize(nb_nodes);
delete dof_synchronizer;
dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension);
dof_synchronizer->initLocalDOFEquationNumbers();
dof_synchronizer->initGlobalDOFEquationNumbers();
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
(*mat_it)->onNodesAdded(nodes_list, event);
}
if (method != _explicit_lumped_mass) {
this->initSolver();
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onNodesRemoved(__attribute__((unused)) const Array<UInt> & element_list,
const Array<UInt> & new_numbering,
__attribute__((unused)) const RemovedNodesEvent & event) {
if(displacement) mesh.removeNodesFromArray(*displacement, new_numbering);
if(mass ) mesh.removeNodesFromArray(*mass , new_numbering);
if(velocity ) mesh.removeNodesFromArray(*velocity , new_numbering);
if(acceleration) mesh.removeNodesFromArray(*acceleration, new_numbering);
if(force ) mesh.removeNodesFromArray(*force , new_numbering);
if(residual ) mesh.removeNodesFromArray(*residual , new_numbering);
if(blocked_dofs) mesh.removeNodesFromArray(*blocked_dofs, new_numbering);
if(increment_acceleration) mesh.removeNodesFromArray(*increment_acceleration, new_numbering);
if(increment) mesh.removeNodesFromArray(*increment , new_numbering);
delete dof_synchronizer;
dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension);
dof_synchronizer->initLocalDOFEquationNumbers();
dof_synchronizer->initGlobalDOFEquationNumbers();
if (method != _explicit_lumped_mass) {
this->initSolver();
}
}
/* -------------------------------------------------------------------------- */
bool SolidMechanicsModel::isInternal(const std::string & field_name,
const ElementKind & element_kind){
bool is_internal = false;
/// check if at least one material contains field_id as an internal
for (UInt m = 0; m < materials.size() && !is_internal; ++m) {
is_internal |= materials[m]->isInternal(field_name, element_kind);
}
return is_internal;
}
/* -------------------------------------------------------------------------- */
ElementTypeMap<UInt> SolidMechanicsModel::getInternalDataPerElem(const std::string & field_name,
const ElementKind & element_kind){
if (!(this->isInternal(field_name,element_kind))) AKANTU_EXCEPTION("unknown internal " << field_name);
for (UInt m = 0; m < materials.size() ; ++m) {
if (materials[m]->isInternal(field_name, element_kind))
return materials[m]->getInternalDataPerElem(field_name,element_kind);
}
return ElementTypeMap<UInt>();
}
/* -------------------------------------------------------------------------- */
ElementTypeMapArray<Real> & SolidMechanicsModel::flattenInternal(const std::string & field_name,
const ElementKind & kind, const GhostType ghost_type){
std::pair<std::string,ElementKind> key(field_name,kind);
if (this->registered_internals.count(key) == 0){
this->registered_internals[key] =
new ElementTypeMapArray<Real>(field_name,this->id);
}
ElementTypeMapArray<Real> * internal_flat = this->registered_internals[key];
for (UInt m = 0; m < materials.size(); ++m)
materials[m]->flattenInternal(field_name,*internal_flat,ghost_type,kind);
return *internal_flat;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::flattenAllRegisteredInternals(const ElementKind & kind){
std::map<std::pair<std::string,ElementKind>,ElementTypeMapArray<Real> *> ::iterator it = this->registered_internals.begin();
std::map<std::pair<std::string,ElementKind>,ElementTypeMapArray<Real> *>::iterator end = this->registered_internals.end();
while (it != end){
if (kind == it->first.second)
this->flattenInternal(it->first.first,kind);
++it;
}
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onDump(){
this->flattenAllRegisteredInternals(_ek_regular);
}
/* -------------------------------------------------------------------------- */
#ifdef AKANTU_USE_IOHELPER
dumper::Field * SolidMechanicsModel
::createElementalField(const std::string & field_name,
const std::string & group_name,
bool padding_flag,
const ElementKind & kind){
dumper::Field * field = NULL;
if(field_name == "partitions")
field = mesh.createElementalField<UInt, dumper::ElementPartitionField>(mesh.getConnectivities(),group_name,this->spatial_dimension,kind);
else if(field_name == "material_index")
field = mesh.createElementalField<UInt, Vector, dumper::ElementalField >(material_index,group_name,this->spatial_dimension,kind);
else {
+ // this copy of field_name is used to compute derivated data such as
+ // strain and von mises stress that are based on grad_u and stress
+ std::string field_name_copy(field_name);
- bool is_internal = this->isInternal(field_name,kind);
+ if (field_name == "strain")
+ field_name_copy = "grad_u";
+ else if (field_name == "Von Mises stress")
+ field_name_copy = "stress";
+
+ bool is_internal = this->isInternal(field_name_copy,kind);
if (is_internal) {
- ElementTypeMap<UInt> nb_data_per_elem = this->getInternalDataPerElem(field_name,kind);
- ElementTypeMapArray<Real> & internal_flat = this->flattenInternal(field_name,kind);
+ ElementTypeMap<UInt> nb_data_per_elem = this->getInternalDataPerElem(field_name_copy,kind);
+ ElementTypeMapArray<Real> & internal_flat = this->flattenInternal(field_name_copy,kind);
field = mesh.createElementalField<Real, dumper::InternalMaterialField>(internal_flat,
group_name,
this->spatial_dimension,kind,nb_data_per_elem);
+ if (field_name == "strain"){
+ dumper::ComputeStrain * foo = new dumper::ComputeStrain(*this);
+ field = dumper::FieldComputeProxy::createFieldCompute(field,*foo);
+ } else if (field_name == "Von Mises stress") {
+ dumper::ComputeVonMisesStress * foo = new dumper::ComputeVonMisesStress(*this);
+ field = dumper::FieldComputeProxy::createFieldCompute(field,*foo);
+ }
//treat the paddings
if (padding_flag){
if (field_name == "stress"){
if (this->spatial_dimension == 2) {
dumper::StressPadder<2> * foo = new dumper::StressPadder<2>(*this);
field = dumper::FieldComputeProxy::createFieldCompute(field,*foo);
}
+ } else if (field_name == "strain"){
+ if (this->spatial_dimension == 2) {
+ dumper::StrainPadder<2> * foo = new dumper::StrainPadder<2>(*this);
+ field = dumper::FieldComputeProxy::createFieldCompute(field,*foo);
+ }
}
- // else if (field_name == "strain"){
- // if (this->spatial_dimension == 2) {
- // dumper::StrainPadder<2> * foo = new dumper::StrainPadder<2>(*this);
- // field = dumper::FieldComputeProxy::createFieldCompute(field,*foo);
- // }
- // }
}
// homogenize the field
dumper::ComputeFunctorInterface * foo =
dumper::HomogenizerProxy::createHomogenizer(*field);
field = dumper::FieldComputeProxy::createFieldCompute(field,*foo);
}
}
return field;
}
/* -------------------------------------------------------------------------- */
dumper::Field * SolidMechanicsModel::createNodalFieldReal(const std::string & field_name,
const std::string & group_name,
bool padding_flag) {
std::map<std::string,Array<Real>* > real_nodal_fields;
real_nodal_fields["displacement" ] = displacement;
real_nodal_fields["mass" ] = mass;
real_nodal_fields["velocity" ] = velocity;
real_nodal_fields["acceleration" ] = acceleration;
real_nodal_fields["force" ] = force;
real_nodal_fields["residual" ] = residual;
real_nodal_fields["increment" ] = increment;
dumper::Field * field = NULL;
if (padding_flag)
field = mesh.createNodalField(real_nodal_fields[field_name],group_name, 3);
else
field = mesh.createNodalField(real_nodal_fields[field_name],group_name);
return field;
}
/* -------------------------------------------------------------------------- */
dumper::Field * SolidMechanicsModel::createNodalFieldBool(const std::string & field_name,
const std::string & group_name,
bool padding_flag) {
std::map<std::string,Array<bool>* > uint_nodal_fields;
uint_nodal_fields["blocked_dofs" ] = blocked_dofs;
dumper::Field * field = NULL;
field = mesh.createNodalField(uint_nodal_fields[field_name],group_name);
return field;
}
/* -------------------------------------------------------------------------- */
#else
/* -------------------------------------------------------------------------- */
dumper::Field * SolidMechanicsModel
::createElementalField(const std::string & field_name,
const std::string & group_name,
bool padding_flag,
const ElementKind & kind){
return NULL;
}
/* -------------------------------------------------------------------------- */
dumper::Field * SolidMechanicsModel::createNodalFieldReal(const std::string & field_name,
const std::string & group_name,
bool padding_flag) {
return NULL;
}
/* -------------------------------------------------------------------------- */
dumper::Field * SolidMechanicsModel::createNodalFieldBool(const std::string & field_name,
const std::string & group_name,
bool padding_flag) {
return NULL;
}
#endif
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::dump(const std::string & dumper_name) {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
synch_registry->synchronize(_gst_for_dump);
mesh.dump(dumper_name);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::dump(const std::string & dumper_name, UInt step) {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
synch_registry->synchronize(_gst_for_dump);
mesh.dump(dumper_name, step);
}
/* ------------------------------------------------------------------------- */
void SolidMechanicsModel::dump(const std::string & dumper_name, Real time, UInt step) {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
synch_registry->synchronize(_gst_for_dump);
mesh.dump(dumper_name, time, step);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::dump() {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
mesh.dump();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::dump(UInt step) {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
mesh.dump(step);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::dump(Real time, UInt step) {
this->onDump();
EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent());
mesh.dump(time, step);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::computeCauchyStresses() {
AKANTU_DEBUG_IN();
// call compute stiffness matrix on each local elements
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
Material & mat = **mat_it;
if(mat.isFiniteDeformation())
mat.computeAllCauchyStresses(_not_ghost);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::saveStressAndStrainBeforeDamage() {
EventManager::sendEvent(SolidMechanicsModelEvent::BeginningOfDamageIterationEvent());
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::updateEnergiesAfterDamage() {
EventManager::sendEvent(SolidMechanicsModelEvent::AfterDamageEvent());
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::printself(std::ostream & stream, int indent) const {
std::string space;
for(Int i = 0; i < indent; i++, space += AKANTU_INDENT);
stream << space << "Solid Mechanics Model [" << std::endl;
stream << space << " + id : " << id << std::endl;
stream << space << " + spatial dimension : " << spatial_dimension << std::endl;
stream << space << " + fem [" << std::endl;
getFEEngine().printself(stream, indent + 2);
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << " + nodals information [" << std::endl;
displacement->printself(stream, indent + 2);
mass ->printself(stream, indent + 2);
velocity ->printself(stream, indent + 2);
acceleration->printself(stream, indent + 2);
force ->printself(stream, indent + 2);
residual ->printself(stream, indent + 2);
blocked_dofs->printself(stream, indent + 2);
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << " + material information [" << std::endl;
material_index.printself(stream, indent + 2);
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << " + materials [" << std::endl;
std::vector<Material *>::const_iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
const Material & mat = *(*mat_it);
mat.printself(stream, indent + 1);
}
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
__END_AKANTU__