diff --git a/packages/20_cohesive_element.cmake b/packages/20_cohesive_element.cmake
index 486c6ce0d..eabcf9e81 100644
--- a/packages/20_cohesive_element.cmake
+++ b/packages/20_cohesive_element.cmake
@@ -1,88 +1,90 @@
#===============================================================================
# @file cohesive_element.cmake
#
# @author Nicolas Richart <nicolas.richart@epfl.ch>
#
# @date Tue Oct 16 14:05:02 2012
#
# @brief package description for cohesive elements
#
# @section LICENSE
#
# Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# Akantu is free software: you can redistribute it and/or modify it under the
# terms of the GNU Lesser General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
# A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
# details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with Akantu. If not, see <http://www.gnu.org/licenses/>.
#
#===============================================================================
option(AKANTU_COHESIVE_ELEMENT "Use cohesive_element package of Akantu" OFF)
add_external_package_dependencies(cohesive_element lapack)
set(AKANTU_COHESIVE_ELEMENT_FILES
model/solid_mechanics/materials/material_cohesive_includes.hh
mesh_utils/cohesive_element_inserter.hh
mesh_utils/cohesive_element_inserter.cc
fem/cohesive_element.cc
fem/shape_cohesive.hh
fem/cohesive_element.hh
fem/fem_template_cohesive.cc
fem/shape_cohesive_inline_impl.cc
model/solid_mechanics/materials/material_cohesive/cohesive_internal_field_tmpl.hh
model/solid_mechanics/materials/material_cohesive/cohesive_internal_field.hh
model/solid_mechanics/materials/material_cohesive/material_cohesive_inline_impl.cc
model/solid_mechanics/solid_mechanics_model_cohesive.cc
+ model/solid_mechanics/fragment_manager.cc
model/solid_mechanics/materials/material_cohesive/material_cohesive.cc
model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.cc
model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_bilinear.cc
model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.cc
model/solid_mechanics/solid_mechanics_model_cohesive.hh
+ model/solid_mechanics/fragment_manager.hh
model/solid_mechanics/materials/material_cohesive/material_cohesive.hh
model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_bilinear.hh
model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.hh
model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.hh
io/dumper/dumper_iohelper_tmpl_connectivity_field.hh
)
set(AKANTU_COHESIVE_ELEMENT_TESTS
test_cohesive_buildfacets_tetrahedron
test_cohesive_buildfacets_hexahedron
test_cohesive_intrinsic
test_cohesive_intrinsic_quadrangle
test_cohesive_extrinsic
test_cohesive_extrinsic_quadrangle
test_cohesive_buildfragments
test_cohesive_intrinsic_impl
)
set(AKANTU_COHESIVE_ELEMENT_DOC
manual/manual-cohesive_element.tex
)
set(AKANTU_COHESIVE_ELEMENT_MANUAL_FILES
manual-cohesive_element.tex
figures/cohesive2d.pdf
figures/cohesive3d.pdf
figures/cohesive_exponential.pdf
)
diff --git a/src/fem/element_group.hh b/src/fem/element_group.hh
index 5ab6271c0..2198ff7a2 100644
--- a/src/fem/element_group.hh
+++ b/src/fem/element_group.hh
@@ -1,153 +1,153 @@
/**
* @file element_group.hh
*
* @author Dana Christen <dana.christen@gmail.com>
*
* @date Wed Mar 06 09:30:00 2013
*
* @brief Stores information relevent to the notion of domain boundary and surfaces.
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_ELEMENT_GROUP_HH__
#define __AKANTU_ELEMENT_GROUP_HH__
#include <set>
#include "aka_common.hh"
#include "aka_memory.hh"
#include "by_element_type.hh"
#include "node_group.hh"
#include "dumpable.hh"
__BEGIN_AKANTU__
class Mesh;
class Element;
/* -------------------------------------------------------------------------- */
class ElementGroup : private Memory, public Dumpable {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
ElementGroup(const std::string & name,
const Mesh & mesh,
NodeGroup & node_group,
UInt dimension = _all_dimensions,
const std::string & id = "element_group",
const MemoryID & memory_id = 0);
/* ------------------------------------------------------------------------ */
/* Type definitions */
/* ------------------------------------------------------------------------ */
public:
typedef ByElementTypeArray<UInt> ElementList;
typedef Array<UInt> NodeList;
/* ------------------------------------------------------------------------ */
/* Node iterator */
/* ------------------------------------------------------------------------ */
typedef NodeGroup::const_node_iterator const_node_iterator;
inline const_node_iterator node_begin() const;
inline const_node_iterator node_end() const;
/* ------------------------------------------------------------------------ */
/* Element iterator */
/* ------------------------------------------------------------------------ */
typedef ElementList::type_iterator type_iterator;
inline type_iterator firstType(UInt dim = _all_dimensions,
const GhostType & ghost_type = _not_ghost,
const ElementKind & kind = _ek_regular) const;
inline type_iterator lastType(UInt dim = _all_dimensions,
const GhostType & ghost_type = _not_ghost,
const ElementKind & kind = _ek_regular) const;
typedef Array<UInt>::const_iterator<UInt> const_element_iterator;
inline const_element_iterator element_begin(const ElementType & type,
- const GhostType & ghost_type) const;
+ const GhostType & ghost_type = _not_ghost) const;
inline const_element_iterator element_end(const ElementType & type,
- const GhostType & ghost_type) const;
+ const GhostType & ghost_type = _not_ghost) const;
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// empty the element group
void empty();
/// append another group to this group
/// BE CAREFUL: it doesn't conserve the element order
void append(const ElementGroup & other_group);
inline void add(const Element & el);
inline void addNode(UInt node_id);
/// function to print the contain of the class
virtual void printself(std::ostream & stream, int indent = 0) const;
private:
inline void addElement(const ElementType & elem_type,
UInt elem_id,
const GhostType & ghost_type);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Elements, elements, UInt);
AKANTU_GET_MACRO(Elements, elements, const ByElementTypeArray<UInt> &);
AKANTU_GET_MACRO(Nodes, node_group.getNodes(), const Array<UInt> &);
AKANTU_GET_MACRO(NodeGroup, node_group, const NodeGroup &);
AKANTU_GET_MACRO_NOT_CONST(NodeGroup, node_group, NodeGroup &);
AKANTU_GET_MACRO(Dimension, dimension, UInt);
AKANTU_GET_MACRO(Name, name, std::string);
inline UInt getNbNodes() const;
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
/// Mesh to which this group belongs
const Mesh & mesh;
/// name of the group
std::string name;
/// list of elements composing the group
ElementList elements;
/// sub list of nodes which are composing the elements
NodeGroup & node_group;
/// group dimension
UInt dimension;
};
#include "element_group_inline_impl.cc"
__END_AKANTU__
#endif /* __AKANTU_ELEMENT_GROUP_HH__ */
diff --git a/src/fem/fem.hh b/src/fem/fem.hh
index 35fbe3ca6..be56d53a8 100644
--- a/src/fem/fem.hh
+++ b/src/fem/fem.hh
@@ -1,360 +1,363 @@
/**
* @file fem.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Tue Jul 20 23:40:43 2010
*
* @brief FEM class
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_FEM_HH__
#define __AKANTU_FEM_HH__
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "aka_memory.hh"
#include "mesh.hh"
#include "element_class.hh"
#include "sparse_matrix.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
class Integrator;
class ShapeFunctions;
}
__BEGIN_AKANTU__
class QuadraturePoint : public Element {
public:
QuadraturePoint(ElementType type = _not_defined, UInt element = 0,
UInt num_point = 0, GhostType ghost_type = _not_ghost) :
Element(type, element, ghost_type), num_point(num_point), global_num(0),
position((Real *)NULL, 0) { };
QuadraturePoint(UInt element, UInt num_point,
UInt global_num,
const Vector<Real> & position,
ElementType type,
GhostType ghost_type = _not_ghost) :
Element(type, element, ghost_type), num_point(num_point), global_num(global_num),
position((Real *)NULL, 0) { this->position.shallowCopy(position); };
QuadraturePoint(const QuadraturePoint & quad) :
Element(quad), num_point(quad.num_point), global_num(quad.global_num), position((Real *) NULL, 0) {
position.shallowCopy(quad.position);
};
inline QuadraturePoint & operator=(const QuadraturePoint & q) {
if(this != &q) {
element = q.element;
type = q.type;
ghost_type = q.ghost_type;
num_point = q.num_point;
global_num = q.global_num;
position.shallowCopy(q.position);
}
return *this;
}
AKANTU_GET_MACRO(Position, position, const Vector<Real> &);
void setPosition(const Vector<Real> & position) {
this->position.shallowCopy(position);
}
/// function to print the containt of the class
virtual void printself(std::ostream & stream, int indent = 0) const {
std::string space;
for(Int i = 0; i < indent; i++, space += AKANTU_INDENT);
stream << space << "QuadraturePoint [";
Element::printself(stream, 0);
stream << ", " << num_point << "]";
}
public:
UInt num_point;
UInt global_num;
private:
Vector<Real> position;
};
/**
* The generic FEM class derived in a FEMTemplate class containing the
* shape functions and the integration method
*/
class FEM : protected Memory {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
FEM(Mesh & mesh, UInt spatial_dimension = _all_dimensions,
ID id = "fem", MemoryID memory_id = 0);
virtual ~FEM();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// build the profile of the sparse matrix corresponding to the mesh
void initSparseMatrixProfile(SparseMatrixType sparse_matrix_type = _unsymmetric);
/// pre-compute all the shape functions, their derivatives and the jacobians
virtual void initShapeFunctions(const GhostType & ghost_type = _not_ghost) = 0;
/// extract the nodal values and store them per element
template<typename T>
static void extractNodalToElementField(const Mesh & mesh,
const Array<T> & nodal_f,
Array<T> & elemental_f,
const ElementType & type,
const GhostType & ghost_type = _not_ghost,
const Array<UInt> & filter_elements = empty_filter);
/// filter a field
template<typename T>
static void filterElementalData(const Mesh & mesh,
const Array<T> & quad_f,
Array<T> & filtered_f,
const ElementType & type,
const GhostType & ghost_type = _not_ghost,
const Array<UInt> & filter_elements = empty_filter);
/* ------------------------------------------------------------------------ */
/* Integration method bridges */
/* ------------------------------------------------------------------------ */
/// integrate f for all elements of type "type"
virtual void integrate(const Array<Real> & f,
Array<Real> &intf,
UInt nb_degree_of_freedom,
const ElementType & type,
const GhostType & ghost_type = _not_ghost,
const Array<UInt> & filter_elements = empty_filter) const = 0;
/// integrate a scalar value on all elements of type "type"
virtual Real integrate(const Array<Real> & f,
const ElementType & type,
const GhostType & ghost_type = _not_ghost,
const Array<UInt> & filter_elements = empty_filter) const = 0;
/// integrate f for all quadrature points of type "type"
virtual void integrateOnQuadraturePoints(const Array<Real> & f,
Array<Real> &intf,
UInt nb_degree_of_freedom,
const ElementType & type,
const GhostType & ghost_type = _not_ghost,
const Array<UInt> & filter_elements = empty_filter) const = 0;
/// integrate one element scalar value on all elements of type "type"
virtual Real integrate(const Vector<Real> & f,
const ElementType & type,
UInt index, const GhostType & ghost_type = _not_ghost) const = 0;
/* ------------------------------------------------------------------------ */
/* compatibility with old FEM fashion */
/* ------------------------------------------------------------------------ */
/// get the number of quadrature points
virtual UInt getNbQuadraturePoints(const ElementType & type,
const GhostType & ghost_type = _not_ghost) const = 0;
/// get the precomputed shapes
const virtual Array<Real> & getShapes(const ElementType & type,
const GhostType & ghost_type = _not_ghost) const = 0;
/// get the derivatives of shapes
const virtual Array<Real> & getShapesDerivatives(const ElementType & type,
const GhostType & ghost_type = _not_ghost,
UInt id = 0) const = 0;
/// get quadrature points
const virtual Matrix<Real> & getQuadraturePoints(const ElementType & type,
const GhostType & ghost_type = _not_ghost) const = 0;
/* ------------------------------------------------------------------------ */
/* Shape method bridges */
/* ------------------------------------------------------------------------ */
virtual
void gradientOnQuadraturePoints(const Array<Real> &u,
Array<Real> &nablauq,
const UInt nb_degree_of_freedom,
const ElementType & type,
const GhostType & ghost_type = _not_ghost,
const Array<UInt> & filter_elements = empty_filter) const = 0;
virtual
void interpolateOnQuadraturePoints(const Array<Real> &u,
Array<Real> &uq,
UInt nb_degree_of_freedom,
const ElementType & type,
const GhostType & ghost_type = _not_ghost,
- const Array<UInt> & filter_elements = empty_filter) const =0;
-
+ const Array<UInt> & filter_elements = empty_filter) const = 0;
+ virtual
+ void interpolateOnQuadraturePoints(const Array<Real> & u,
+ ByElementTypeReal & uq,
+ const ByElementTypeUInt * filter_elements = NULL) const = 0;
/* ------------------------------------------------------------------------ */
/* Other methods */
/* ------------------------------------------------------------------------ */
/// pre-compute normals on control points
virtual void computeNormalsOnControlPoints(const GhostType & ghost_type = _not_ghost) = 0;
/// pre-compute normals on control points
virtual void computeNormalsOnControlPoints(__attribute__((unused)) const Array<Real> & field,
__attribute__((unused)) const GhostType & ghost_type = _not_ghost) {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/// pre-compute normals on control points
virtual void computeNormalsOnControlPoints(__attribute__((unused)) const Array<Real> & field,
__attribute__((unused)) Array<Real> & normal,
__attribute__((unused)) const ElementType & type,
__attribute__((unused)) const GhostType & ghost_type = _not_ghost) const {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/// assemble vectors
void assembleArray(const Array<Real> & elementary_vect,
Array<Real> & nodal_values,
const Array<Int> & equation_number,
UInt nb_degree_of_freedom,
const ElementType & type,
const GhostType & ghost_type = _not_ghost,
const Array<UInt> & filter_elements = empty_filter,
Real scale_factor = 1) const;
/// assemble matrix in the complete sparse matrix
void assembleMatrix(const Array<Real> & elementary_mat,
SparseMatrix & matrix,
UInt nb_degree_of_freedom,
const ElementType & type,
const GhostType & ghost_type = _not_ghost,
const Array<UInt> & filter_elements = empty_filter) const;
/// assemble a field as a lumped matrix (ex. rho in lumped mass)
virtual void assembleFieldLumped(__attribute__ ((unused)) const Array<Real> & field_1,
__attribute__ ((unused)) UInt nb_degree_of_freedom,
__attribute__ ((unused)) Array<Real> & lumped,
__attribute__ ((unused)) const Array<Int> & equation_number,
__attribute__ ((unused)) ElementType type,
__attribute__ ((unused)) const GhostType & ghost_type) const {
AKANTU_DEBUG_TO_IMPLEMENT();
};
/// assemble a field as a matrix (ex. rho to mass matrix)
virtual void assembleFieldMatrix(__attribute__ ((unused)) const Array<Real> & field_1,
__attribute__ ((unused)) UInt nb_degree_of_freedom,
__attribute__ ((unused)) SparseMatrix & matrix,
__attribute__ ((unused)) ElementType type,
__attribute__ ((unused)) const GhostType & ghost_type) const {
AKANTU_DEBUG_TO_IMPLEMENT();
}
/// function to print the containt of the class
virtual void printself(std::ostream & stream, int indent = 0) const;
private:
/// initialise the class
void init();
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
AKANTU_GET_MACRO(ElementDimension, element_dimension, UInt);
/// get the mesh contained in the fem object
inline Mesh & getMesh() const;
/// get the in-radius of an element
static inline Real getElementInradius(const Matrix<Real> & coord, const ElementType & type);
/// get the normals on quadrature points
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(NormalsOnQuadPoints, normals_on_quad_points, Real);
/// get cohesive element type for a given facet type
static inline ElementType getCohesiveElementType(const ElementType & type_facet);
virtual const ShapeFunctions & getShapeFunctionsInterface() const = 0;
virtual const Integrator & getIntegratorInterface() const = 0;
static inline InterpolationType getInterpolationType(const ElementType & el_type);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
/// spatial dimension of the problem
UInt element_dimension;
/// the mesh on which all computation are made
Mesh & mesh;
/// normals at quadrature points
ByElementTypeReal normals_on_quad_points;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#if defined (AKANTU_INCLUDE_INLINE_IMPL)
# include "fem_inline_impl.cc"
#endif
/// standard output stream operator
inline std::ostream & operator <<(std::ostream & stream, const FEM & _this)
{
_this.printself(stream);
return stream;
}
/// standard output stream operator
inline std::ostream & operator <<(std::ostream & stream, const QuadraturePoint & _this)
{
_this.printself(stream);
return stream;
}
__END_AKANTU__
#include "fem_template.hh"
#endif /* __AKANTU_FEM_HH__ */
diff --git a/src/fem/fem_template.hh b/src/fem/fem_template.hh
index 0e64883b3..37ec8647f 100644
--- a/src/fem/fem_template.hh
+++ b/src/fem/fem_template.hh
@@ -1,255 +1,260 @@
/**
* @file fem_template.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Tue Feb 15 16:32:44 2011
*
* @brief templated class that calls integration and shape objects
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
#ifndef __AKANTU_FEM_TEMPLATE_HH__
#define __AKANTU_FEM_TEMPLATE_HH__
/* -------------------------------------------------------------------------- */
#include "fem.hh"
#include "integrator.hh"
#include "shape_functions.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind = _ek_regular>
class FEMTemplate : public FEM {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
typedef I<kind> Integ;
typedef S<kind> Shape;
FEMTemplate(Mesh & mesh, UInt spatial_dimension = _all_dimensions,
ID id = "fem", MemoryID memory_id = 0);
virtual ~FEMTemplate();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// pre-compute all the shape functions, their derivatives and the jacobians
void initShapeFunctions(const GhostType & ghost_type = _not_ghost);
void initShapeFunctions(const Array<Real> & nodes,
const GhostType & ghost_type = _not_ghost);
/* ------------------------------------------------------------------------ */
/* Integration method bridges */
/* ------------------------------------------------------------------------ */
/// integrate f for all elements of type "type"
void integrate(const Array<Real> & f,
Array<Real> &intf,
UInt nb_degree_of_freedom,
const ElementType & type,
const GhostType & ghost_type = _not_ghost,
const Array<UInt> & filter_elements = empty_filter) const;
/// integrate a scalar value on all elements of type "type"
Real integrate(const Array<Real> & f,
const ElementType & type,
const GhostType & ghost_type = _not_ghost,
const Array<UInt> & filter_elements = empty_filter) const;
/// integrate one element scalar value on all elements of type "type"
virtual Real integrate(const Vector<Real> & f,
const ElementType & type,
UInt index, const GhostType & ghost_type = _not_ghost) const;
/// integrate partially around a quadrature point (@f$ intf_q = f_q * J_q * w_q @f$)
void integrateOnQuadraturePoints(const Array<Real> & f,
Array<Real> &intf,
UInt nb_degree_of_freedom,
const ElementType & type,
const GhostType & ghost_type = _not_ghost,
const Array<UInt> & filter_elements = empty_filter) const;
/// get the number of quadrature points
UInt getNbQuadraturePoints(const ElementType & type,
const GhostType & ghost_type = _not_ghost) const;
/// get shapes precomputed
const Array<Real> & getShapes(const ElementType & type,
- const GhostType & ghost_type = _not_ghost) const;
+ const GhostType & ghost_type = _not_ghost) const;
/// get the derivatives of shapes
const Array<Real> & getShapesDerivatives(const ElementType & type,
- const GhostType & ghost_type = _not_ghost,
- UInt id=0) const;
+ const GhostType & ghost_type = _not_ghost,
+ UInt id=0) const;
/// get quadrature points
const inline Matrix<Real> & getQuadraturePoints(const ElementType & type,
- const GhostType & ghost_type = _not_ghost) const;
+ const GhostType & ghost_type = _not_ghost) const;
/* ------------------------------------------------------------------------ */
/* Shape method bridges */
/* ------------------------------------------------------------------------ */
/// compute the gradient of a nodal field on the quadrature points
void gradientOnQuadraturePoints(const Array<Real> &u,
Array<Real> &nablauq,
const UInt nb_degree_of_freedom,
const ElementType & type,
const GhostType & ghost_type = _not_ghost,
const Array<UInt> & filter_elements = empty_filter) const;
/// interpolate a nodal field on the quadrature points
void interpolateOnQuadraturePoints(const Array<Real> &u,
Array<Real> &uq,
UInt nb_degree_of_freedom,
const ElementType & type,
const GhostType & ghost_type = _not_ghost,
const Array<UInt> & filter_elements = empty_filter) const;
+ /// interpolate a nodal field on the quadrature points given a by_element_type
+ void interpolateOnQuadraturePoints(const Array<Real> & u,
+ ByElementTypeReal & uq,
+ const ByElementTypeUInt * filter_elements = NULL) const;
+
/// find natural coords from real coords provided an element
void inverseMap(const Vector<Real> & real_coords,
UInt element,
const ElementType & type,
Vector<Real> & natural_coords,
const GhostType & ghost_type = _not_ghost) const;
/// return true if the coordinates provided are inside the element, false otherwise
inline bool contains(const Vector<Real> & real_coords,
UInt element,
const ElementType & type,
const GhostType & ghost_type = _not_ghost) const;
/// compute the shape on a provided point
inline void computeShapes(const Vector<Real> & real_coords,
UInt element,
const ElementType & type,
Vector<Real> & shapes,
const GhostType & ghost_type = _not_ghost) const;
/* ------------------------------------------------------------------------ */
/* Other methods */
/* ------------------------------------------------------------------------ */
/// pre-compute normals on control points
void computeNormalsOnControlPoints(const GhostType & ghost_type = _not_ghost);
void computeNormalsOnControlPoints(const Array<Real> & field,
const GhostType & ghost_type = _not_ghost);
void computeNormalsOnControlPoints(const Array<Real> & field,
Array<Real> & normal,
const ElementType & type,
const GhostType & ghost_type = _not_ghost) const;
template<ElementType type>
void computeNormalsOnControlPoints(const Array<Real> & field,
Array<Real> & normal,
const GhostType & ghost_type) const;
/// function to print the contain of the class
// virtual void printself(std::ostream & stream, int indent = 0) const{};
void assembleFieldLumped(const Array<Real> & field_1,
UInt nb_degree_of_freedom,
Array<Real> & lumped,
const Array<Int> & equation_number,
ElementType type,
const GhostType & ghost_type = _not_ghost) const;
void assembleFieldMatrix(const Array<Real> & field,
UInt nb_degree_of_freedom,
SparseMatrix & matrix,
ElementType type,
const GhostType & ghost_type = _not_ghost) const;
private:
template <ElementType type>
void assembleLumpedTemplate(const Array<Real> & field_1,
UInt nb_degree_of_freedom,
Array<Real> & lumped,
const Array<Int> & equation_number,
const GhostType & ghost_type) const;
/// @f$ \tilde{M}_{i} = \sum_j M_{ij} = \sum_j \int \rho \varphi_i \varphi_j dV = \int \rho \varphi_i dV @f$
template <ElementType type>
void assembleLumpedRowSum(const Array<Real> & field_1,
UInt nb_degree_of_freedom,
Array<Real> & lumped,
const Array<Int> & equation_number,
const GhostType & ghost_type) const;
/// @f$ \tilde{M}_{i} = c * M_{ii} = \int_{V_e} \rho dV @f$
template <ElementType type>
void assembleLumpedDiagonalScaling(const Array<Real> & field_1,
UInt nb_degree_of_freedom,
Array<Real> & lumped,
const Array<Int> & equation_number,
const GhostType & ghost_type) const;
template <ElementType type>
void assembleFieldMatrix(const Array<Real> & field,
UInt nb_degree_of_freedom,
SparseMatrix & matrix,
const GhostType & ghost_type) const;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
const ShapeFunctions & getShapeFunctionsInterface() const { return shape_functions; };
const Shape & getShapeFunctions() const { return shape_functions; };
const Integrator & getIntegratorInterface() const { return integrator; };
const Integ & getIntegrator() const { return integrator; };
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
Integ integrator;
Shape shape_functions;
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#include "fem_template_tmpl.hh"
__END_AKANTU__
#endif /* __AKANTU_FEM_TEMPLATE_HH__ */
diff --git a/src/fem/fem_template_tmpl.hh b/src/fem/fem_template_tmpl.hh
index e408df674..f6c7399aa 100644
--- a/src/fem/fem_template_tmpl.hh
+++ b/src/fem/fem_template_tmpl.hh
@@ -1,983 +1,1034 @@
/**
* @file fem_template_tmpl.hh
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @date Mon Nov 5 17:08:50 2012
*
* @brief Template implementation of FEMTemplate
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
FEMTemplate<I, S, kind>::FEMTemplate(Mesh & mesh, UInt spatial_dimension,
ID id, MemoryID memory_id) :
FEM(mesh,spatial_dimension,id,memory_id),
integrator(mesh, id, memory_id),
shape_functions(mesh, id, memory_id) { }
/* -------------------------------------------------------------------------- */
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
FEMTemplate<I, S, kind>::~FEMTemplate() { }
/* -------------------------------------------------------------------------- */
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
void FEMTemplate<I, S, kind>::gradientOnQuadraturePoints(const Array<Real> &u,
Array<Real> &nablauq,
const UInt nb_degree_of_freedom,
const ElementType & type,
const GhostType & ghost_type,
const Array<UInt> & filter_elements) const {
AKANTU_DEBUG_IN();
UInt nb_element = mesh.getNbElement(type, ghost_type);
if(filter_elements != empty_filter) nb_element = filter_elements.getSize();
UInt nb_points = shape_functions.getControlPoints(type, ghost_type).cols();
#ifndef AKANTU_NDEBUG
UInt element_dimension = mesh.getSpatialDimension(type);
AKANTU_DEBUG_ASSERT(u.getSize() == mesh.getNbNodes(),
"The vector u(" << u.getID()
<< ") has not the good size.");
AKANTU_DEBUG_ASSERT(u.getNbComponent() == nb_degree_of_freedom ,
"The vector u(" << u.getID()
<< ") has not the good number of component.");
AKANTU_DEBUG_ASSERT(nablauq.getNbComponent()
== nb_degree_of_freedom * element_dimension,
"The vector nablauq(" << nablauq.getID()
<< ") has not the good number of component.");
// AKANTU_DEBUG_ASSERT(nablauq.getSize() == nb_element * nb_points,
- // "The vector nablauq(" << nablauq.getID()
- // << ") has not the good size.");
+ // "The vector nablauq(" << nablauq.getID()
+ // << ") has not the good size.");
#endif
nablauq.resize(nb_element * nb_points);
#define COMPUTE_GRADIENT(type) \
if (element_dimension == ElementClass<type>::getSpatialDimension()) \
shape_functions.template gradientOnControlPoints<type>(u, \
nablauq, \
nb_degree_of_freedom, \
ghost_type, \
filter_elements);
if(kind == _ek_regular)
AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(COMPUTE_GRADIENT);
else if(kind == _ek_cohesive)
AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_GRADIENT);
else
AKANTU_DEBUG_TO_IMPLEMENT();
#undef COMPUTE_GRADIENT
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
void FEMTemplate<I, S, kind>::initShapeFunctions(const GhostType & ghost_type) {
initShapeFunctions(mesh.getNodes(), ghost_type);
}
/* -------------------------------------------------------------------------- */
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
void FEMTemplate<I, S, kind>::initShapeFunctions(const Array<Real> & nodes,
const GhostType & ghost_type) {
AKANTU_DEBUG_IN();
Mesh::type_iterator it = mesh.firstType(element_dimension, ghost_type, kind);
Mesh::type_iterator end = mesh.lastType(element_dimension, ghost_type, kind);
for(; it != end; ++it) {
ElementType type = *it;
integrator.initIntegrator(nodes, type, ghost_type);
const Matrix<Real> control_points =
getQuadraturePoints(type, ghost_type);
shape_functions.initShapeFunctions(nodes, control_points, type, ghost_type);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
void FEMTemplate<I, S, kind>::integrate(const Array<Real> & f,
Array<Real> &intf,
UInt nb_degree_of_freedom,
const ElementType & type,
const GhostType & ghost_type,
const Array<UInt> & filter_elements) const{
UInt nb_element = mesh.getNbElement(type, ghost_type);
if(filter_elements != empty_filter) nb_element = filter_elements.getSize();
#ifndef AKANTU_NDEBUG
UInt nb_quadrature_points = getNbQuadraturePoints(type);
AKANTU_DEBUG_ASSERT(f.getSize() == nb_element * nb_quadrature_points,
"The vector f(" << f.getID() << " size " << f.getSize()
<< ") has not the good size (" << nb_element << ").");
AKANTU_DEBUG_ASSERT(f.getNbComponent() == nb_degree_of_freedom ,
"The vector f(" << f.getID()
<< ") has not the good number of component.");
AKANTU_DEBUG_ASSERT(intf.getNbComponent() == nb_degree_of_freedom,
"The vector intf(" << intf.getID()
<< ") has not the good number of component.");
#endif
intf.resize(nb_element);
#define INTEGRATE(type) \
integrator.template integrate<type>(f, \
intf, \
nb_degree_of_freedom, \
ghost_type, \
filter_elements);
if(kind == _ek_regular)
AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(INTEGRATE);
else if(kind == _ek_cohesive)
AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(INTEGRATE);
else if(kind == _ek_structural)
AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(INTEGRATE);
else
AKANTU_DEBUG_TO_IMPLEMENT();
#undef INTEGRATE
}
/* -------------------------------------------------------------------------- */
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
Real FEMTemplate<I, S, kind>::integrate(const Array<Real> & f,
const ElementType & type,
const GhostType & ghost_type,
const Array<UInt> & filter_elements) const{
AKANTU_DEBUG_IN();
#ifndef AKANTU_NDEBUG
// std::stringstream sstr; sstr << ghost_type;
// AKANTU_DEBUG_ASSERT(sstr.str() == nablauq.getTag(),
// "The vector " << nablauq.getID() << " is not taged " << ghost_type);
UInt nb_element = mesh.getNbElement(type, ghost_type);
if(filter_elements != empty_filter) nb_element = filter_elements.getSize();
UInt nb_quadrature_points = getNbQuadraturePoints(type, ghost_type);
AKANTU_DEBUG_ASSERT(f.getSize() == nb_element * nb_quadrature_points,
"The vector f(" << f.getID()
<< ") has not the good size. (" << f.getSize() << "!=" << nb_quadrature_points * nb_element << ")");
AKANTU_DEBUG_ASSERT(f.getNbComponent() == 1,
"The vector f(" << f.getID()
<< ") has not the good number of component.");
#endif
Real integral = 0.;
#define INTEGRATE(type) \
integral = integrator.template integrate<type>(f, \
ghost_type, \
filter_elements);
if(kind == _ek_regular)
AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(INTEGRATE);
else if(kind == _ek_cohesive)
AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(INTEGRATE);
else if(kind == _ek_structural)
AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(INTEGRATE);
else
AKANTU_DEBUG_TO_IMPLEMENT();
#undef INTEGRATE
AKANTU_DEBUG_OUT();
return integral;
}
/* -------------------------------------------------------------------------- */
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
Real FEMTemplate<I, S, kind>::integrate(const Vector<Real> & f,
const ElementType & type,
UInt index,
const GhostType & ghost_type) const{
Real res = 0.;
#define INTEGRATE(type) \
res = integrator.template integrate<type>(f, \
index, \
ghost_type);
if(kind == _ek_regular)
AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(INTEGRATE);
else if(kind == _ek_cohesive)
AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(INTEGRATE);
else if(kind == _ek_structural)
AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(INTEGRATE);
else
AKANTU_DEBUG_TO_IMPLEMENT();
#undef INTEGRATE
return res;
}
/* -------------------------------------------------------------------------- */
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
void FEMTemplate<I, S, kind>::integrateOnQuadraturePoints(const Array<Real> & f,
Array<Real> &intf,
UInt nb_degree_of_freedom,
const ElementType & type,
const GhostType & ghost_type,
const Array<UInt> & filter_elements) const{
UInt nb_element = mesh.getNbElement(type, ghost_type);
if(filter_elements != empty_filter) nb_element = filter_elements.getSize();
UInt nb_quadrature_points = getNbQuadraturePoints(type);
#ifndef AKANTU_NDEBUG
// std::stringstream sstr; sstr << ghost_type;
// AKANTU_DEBUG_ASSERT(sstr.str() == nablauq.getTag(),
// "The vector " << nablauq.getID() << " is not taged " << ghost_type);
AKANTU_DEBUG_ASSERT(f.getSize() == nb_element * nb_quadrature_points,
"The vector f(" << f.getID() << " size " << f.getSize()
<< ") has not the good size (" << nb_element << ").");
AKANTU_DEBUG_ASSERT(f.getNbComponent() == nb_degree_of_freedom ,
"The vector f(" << f.getID()
<< ") has not the good number of component.");
AKANTU_DEBUG_ASSERT(intf.getNbComponent() == nb_degree_of_freedom,
"The vector intf(" << intf.getID()
<< ") has not the good number of component.");
#endif
intf.resize(nb_element*nb_quadrature_points);
#define INTEGRATE(type) \
integrator.template integrateOnQuadraturePoints<type>(f, \
intf, \
nb_degree_of_freedom, \
ghost_type, \
filter_elements);
if(kind == _ek_regular)
AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(INTEGRATE);
else if(kind == _ek_cohesive)
AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(INTEGRATE);
else if(kind == _ek_structural)
AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(INTEGRATE);
else
AKANTU_DEBUG_TO_IMPLEMENT();
#undef INTEGRATE
}
/* -------------------------------------------------------------------------- */
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
void FEMTemplate<I, S, kind>::interpolateOnQuadraturePoints(const Array<Real> &u,
Array<Real> &uq,
UInt nb_degree_of_freedom,
const ElementType & type,
const GhostType & ghost_type,
const Array<UInt> & filter_elements) const{
AKANTU_DEBUG_IN();
UInt nb_points = shape_functions.getControlPoints(type, ghost_type).cols();
UInt nb_element = mesh.getNbElement(type, ghost_type);
if(filter_elements != empty_filter) nb_element = filter_elements.getSize();
#ifndef AKANTU_NDEBUG
// std::stringstream sstr; sstr << ghost_type;
// AKANTU_DEBUG_ASSERT(sstr.str() == nablauq.getTag(),
// "The vector " << nablauq.getID() << " is not taged " << ghost_type);
AKANTU_DEBUG_ASSERT(u.getSize() == mesh.getNbNodes(),
"The vector u(" << u.getID()
<< ") has not the good size.");
AKANTU_DEBUG_ASSERT(u.getNbComponent() == nb_degree_of_freedom ,
"The vector u(" << u.getID()
<< ") has not the good number of component.");
AKANTU_DEBUG_ASSERT(uq.getNbComponent() == nb_degree_of_freedom,
"The vector uq(" << uq.getID()
<< ") has not the good number of component.");
#endif
uq.resize(nb_element * nb_points);
#define INTERPOLATE(type) \
shape_functions.template interpolateOnControlPoints<type>(u, \
uq, \
nb_degree_of_freedom, \
ghost_type, \
filter_elements);
if(kind == _ek_regular)
AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(INTERPOLATE);
else if(kind == _ek_cohesive)
AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(INTERPOLATE);
else
AKANTU_DEBUG_TO_IMPLEMENT();
#undef INTERPOLATE
AKANTU_DEBUG_OUT();
}
+/* -------------------------------------------------------------------------- */
+template<template <ElementKind> class I,
+ template <ElementKind> class S,
+ ElementKind kind>
+void FEMTemplate<I, S, kind>::interpolateOnQuadraturePoints(const Array<Real> & u,
+ ByElementTypeReal & uq,
+ const ByElementTypeUInt * filter_elements) const {
+ AKANTU_DEBUG_IN();
+
+ const Array<UInt> * filter = NULL;
+
+ for (ghost_type_t::iterator gt = ghost_type_t::begin();
+ gt != ghost_type_t::end(); ++gt) {
+
+ GhostType ghost_type = *gt;
+ ByElementTypeReal::type_iterator it = uq.firstType(_all_dimensions, ghost_type, kind);
+ ByElementTypeReal::type_iterator last = uq.lastType(_all_dimensions, ghost_type, kind);
+
+ for (; it != last; ++it) {
+ ElementType type = *it;
+
+ UInt nb_quad_per_element = getNbQuadraturePoints(type, ghost_type);
+
+ UInt nb_element = 0;
+
+ if (filter_elements) {
+ filter = &((*filter_elements)(type, ghost_type));
+ nb_element = filter->getSize();
+ }
+ else {
+ filter = &empty_filter;
+ nb_element = mesh.getNbElement(type, ghost_type);
+ }
+
+ UInt nb_tot_quad = nb_quad_per_element * nb_element;
+
+ Array<Real> & quad = uq(type, ghost_type);
+ quad.resize(nb_tot_quad);
+
+ interpolateOnQuadraturePoints(u,
+ quad,
+ mesh.getSpatialDimension(),
+ type,
+ ghost_type,
+ *filter);
+ }
+ }
+
+ AKANTU_DEBUG_OUT();
+}
+
/* -------------------------------------------------------------------------- */
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
void FEMTemplate<I, S, kind>::computeNormalsOnControlPoints(const GhostType & ghost_type) {
AKANTU_DEBUG_IN();
computeNormalsOnControlPoints(mesh.getNodes(),
ghost_type);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
void FEMTemplate<I, S, kind>::computeNormalsOnControlPoints(const Array<Real> & field,
const GhostType & ghost_type) {
AKANTU_DEBUG_IN();
// Real * coord = mesh.getNodes().values;
UInt spatial_dimension = mesh.getSpatialDimension();
//allocate the normal arrays
Mesh::type_iterator it = mesh.firstType(element_dimension, ghost_type, kind);
Mesh::type_iterator end = mesh.lastType(element_dimension, ghost_type, kind);
for(; it != end; ++it) {
ElementType type = *it;
UInt size = mesh.getNbElement(type, ghost_type);
if(normals_on_quad_points.exists(type, ghost_type)) {
normals_on_quad_points(type, ghost_type).resize(size);
} else {
normals_on_quad_points.alloc(size, spatial_dimension, type, ghost_type);
}
}
//loop over the type to build the normals
it = mesh.firstType(element_dimension, ghost_type, kind);
for(; it != end; ++it) {
Array<Real> & normals_on_quad = normals_on_quad_points(*it, ghost_type);
computeNormalsOnControlPoints(field, normals_on_quad, *it, ghost_type);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
void FEMTemplate<I, S, kind>::computeNormalsOnControlPoints(const Array<Real> & field,
- Array<Real> & normal,
- const ElementType & type,
- const GhostType & ghost_type) const {
+ Array<Real> & normal,
+ const ElementType & type,
+ const GhostType & ghost_type) const {
#define COMPUTE_NORMALS_ON_QUAD(type) \
computeNormalsOnControlPoints<type>(field, normal, ghost_type);
if(kind == _ek_regular)
AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(COMPUTE_NORMALS_ON_QUAD);
else if(kind == _ek_cohesive)
AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_NORMALS_ON_QUAD);
else
AKANTU_DEBUG_TO_IMPLEMENT();
#undef COMPUTE_NORMALS_ON_QUAD
}
/* -------------------------------------------------------------------------- */
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
template<ElementType type>
void FEMTemplate<I, S, kind>::computeNormalsOnControlPoints(const Array<Real> & field,
- Array<Real> & normal,
- const GhostType & ghost_type) const {
+ Array<Real> & normal,
+ const GhostType & ghost_type) const {
AKANTU_DEBUG_IN();
UInt spatial_dimension = mesh.getSpatialDimension();
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
UInt nb_points = getNbQuadraturePoints(type, ghost_type);
UInt nb_element = mesh.getConnectivity(type, ghost_type).getSize();
normal.resize(nb_element * nb_points);
Array<Real>::iterator< Matrix<Real> > normals_on_quad = normal.begin_reinterpret(spatial_dimension,
- nb_points,
- nb_element);
+ nb_points,
+ nb_element);
Array<Real> f_el(0, spatial_dimension * nb_nodes_per_element);
FEM::extractNodalToElementField(mesh, field, f_el, type, ghost_type);
const Matrix<Real> quads =
integrator. template getQuadraturePoints<type>(ghost_type);
Array<Real>::iterator< Matrix<Real> > f_it = f_el.begin(spatial_dimension, nb_nodes_per_element);
for (UInt elem = 0; elem < nb_element; ++elem) {
ElementClass<type>::computeNormalsOnNaturalCoordinates(quads,
*f_it,
*normals_on_quad);
++normals_on_quad;
++f_it;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/* Matrix lumping functions */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
void FEMTemplate<I, S, kind>::assembleFieldLumped(const Array<Real> & field_1,
- UInt nb_degree_of_freedom,
- Array<Real> & lumped,
- const Array<Int> & equation_number,
- ElementType type,
- const GhostType & ghost_type) const {
+ UInt nb_degree_of_freedom,
+ Array<Real> & lumped,
+ const Array<Int> & equation_number,
+ ElementType type,
+ const GhostType & ghost_type) const {
AKANTU_DEBUG_IN();
#define ASSEMBLE_LUMPED(type) \
assembleLumpedTemplate<type>(field_1, nb_degree_of_freedom,lumped, equation_number,ghost_type)
AKANTU_BOOST_ALL_ELEMENT_SWITCH(ASSEMBLE_LUMPED);;
#undef ASSEMBLE_LUMPED
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
void FEMTemplate<I, S, kind>::assembleFieldMatrix(const Array<Real> & field_1,
- UInt nb_degree_of_freedom,
- SparseMatrix & matrix,
- ElementType type,
- const GhostType & ghost_type) const {
+ UInt nb_degree_of_freedom,
+ SparseMatrix & matrix,
+ ElementType type,
+ const GhostType & ghost_type) const {
AKANTU_DEBUG_IN();
#define ASSEMBLE_MATRIX(type) \
assembleFieldMatrix<type>(field_1, nb_degree_of_freedom, \
matrix, \
ghost_type)
AKANTU_BOOST_ALL_ELEMENT_SWITCH(ASSEMBLE_MATRIX);;
#undef ASSEMBLE_MATRIX
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
template <ElementType type>
void FEMTemplate<I, S, kind>::assembleLumpedTemplate(const Array<Real> & field_1,
- UInt nb_degree_of_freedom,
- Array<Real> & lumped,
- const Array<Int> & equation_number,
- const GhostType & ghost_type) const {
+ UInt nb_degree_of_freedom,
+ Array<Real> & lumped,
+ const Array<Int> & equation_number,
+ const GhostType & ghost_type) const {
this->template assembleLumpedRowSum<type>(field_1, nb_degree_of_freedom,lumped, equation_number,ghost_type);
}
/* -------------------------------------------------------------------------- */
/**
* @f$ \tilde{M}_{i} = \sum_j M_{ij} = \sum_j \int \rho \varphi_i \varphi_j dV = \int \rho \varphi_i dV @f$
*/
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
template <ElementType type>
void FEMTemplate<I, S, kind>::assembleLumpedRowSum(const Array<Real> & field_1,
- UInt nb_degree_of_freedom,
- Array<Real> & lumped,
- const Array<Int> & equation_number,
- const GhostType & ghost_type) const {
+ UInt nb_degree_of_freedom,
+ Array<Real> & lumped,
+ const Array<Int> & equation_number,
+ const GhostType & ghost_type) const {
AKANTU_DEBUG_IN();
UInt shapes_size = ElementClass<type>::getShapeSize();
Array<Real> * field_times_shapes = new Array<Real>(0, 1);//shapes_size);
shape_functions.template fieldTimesShapes<type>(field_1, *field_times_shapes, ghost_type);
UInt nb_element = mesh.getNbElement(type, ghost_type);
Array<Real> * int_field_times_shapes = new Array<Real>(nb_element, shapes_size,
- "inte_rho_x_shapes");
+ "inte_rho_x_shapes");
integrator.template integrate<type>(*field_times_shapes, *int_field_times_shapes,
shapes_size, ghost_type, empty_filter);
delete field_times_shapes;
int_field_times_shapes->extendComponentsInterlaced(nb_degree_of_freedom,1);
assembleArray(*int_field_times_shapes, lumped, equation_number,nb_degree_of_freedom, type, ghost_type);
delete int_field_times_shapes;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* @f$ \tilde{M}_{i} = c * M_{ii} = \int_{V_e} \rho dV @f$
*/
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
template <ElementType type>
void FEMTemplate<I, S, kind>::assembleLumpedDiagonalScaling(const Array<Real> & field_1,
- UInt nb_degree_of_freedom,
- Array<Real> & lumped,
- const Array<Int> & equation_number,
- const GhostType & ghost_type) const {
+ UInt nb_degree_of_freedom,
+ Array<Real> & lumped,
+ const Array<Int> & equation_number,
+ const GhostType & ghost_type) const {
AKANTU_DEBUG_IN();
UInt nb_nodes_per_element_p1 = ElementClass<type>::getNbNodesPerElement();
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
UInt nb_quadrature_points = integrator.template getQuadraturePoints<type>(ghost_type).cols();
UInt nb_element = field_1.getSize() / nb_quadrature_points;
Real corner_factor = 0;
Real mid_factor = 0;
if(type == _triangle_6) {
corner_factor = 1./12.;
mid_factor = 1./4.;
}
if (type == _tetrahedron_10) {
corner_factor = 1./32.;
mid_factor = 7./48.;
}
if (type == _quadrangle_8) {
corner_factor = 1./36.;
mid_factor = 8./36.;
}
if (nb_element == 0) {
AKANTU_DEBUG_OUT();
return;
}
/// compute @f$ \int \rho dV = \rho V @f$ for each element
Array<Real> * int_field_1 = new Array<Real>(field_1.getSize(), 1,
"inte_rho_x_1");
integrator.template integrate<type>(field_1, *int_field_1, 1, ghost_type, empty_filter);
/// distribute the mass of the element to the nodes
Array<Real> * lumped_per_node = new Array<Real>(nb_element, nb_nodes_per_element, "mass_per_node");
Real * int_field_1_val = int_field_1->values;
Real * lumped_per_node_val = lumped_per_node->values;
for (UInt e = 0; e < nb_element; ++e) {
Real lmass = *int_field_1_val * corner_factor;
for (UInt n = 0; n < nb_nodes_per_element_p1; ++n)
*lumped_per_node_val++ = lmass; /// corner points
lmass = *int_field_1_val * mid_factor;
for (UInt n = nb_nodes_per_element_p1; n < nb_nodes_per_element; ++n)
*lumped_per_node_val++ = lmass; /// mid points
int_field_1_val++;
}
delete int_field_1;
lumped_per_node->extendComponentsInterlaced(nb_degree_of_freedom,1);
assembleArray(*lumped_per_node, lumped, equation_number, nb_degree_of_freedom, type, ghost_type);
delete lumped_per_node;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* @f$ \tilde{M}_{i} = \sum_j M_{ij} = \sum_j \int \rho \varphi_i \varphi_j dV = \int \rho \varphi_i dV @f$
*/
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
template <ElementType type>
void FEMTemplate<I, S, kind>::assembleFieldMatrix(const Array<Real> & field_1,
- UInt nb_degree_of_freedom,
- SparseMatrix & matrix,
- const GhostType & ghost_type) const {
+ UInt nb_degree_of_freedom,
+ SparseMatrix & matrix,
+ const GhostType & ghost_type) const {
AKANTU_DEBUG_IN();
UInt vect_size = field_1.getSize();
UInt shapes_size = ElementClass<type>::getShapeSize();
UInt lmat_size = nb_degree_of_freedom * shapes_size;
const Array<Real> & shapes = shape_functions.getShapes(type,ghost_type);
Array<Real> * modified_shapes = new Array<Real>(vect_size, lmat_size * nb_degree_of_freedom);
modified_shapes->clear();
Array<Real> * local_mat = new Array<Real>(vect_size, lmat_size * lmat_size);
Array<Real>::iterator< Matrix<Real> > shape_vect = modified_shapes->begin(nb_degree_of_freedom, lmat_size);
Real * sh = shapes.values;
for(UInt q = 0; q < vect_size; ++q) {
Real * msh = shape_vect->storage();
for (UInt d = 0; d < nb_degree_of_freedom; ++d) {
Real * msh_tmp = msh + d * (lmat_size + 1);
for (UInt s = 0; s < shapes_size; ++s) {
*msh_tmp = sh[s];
msh_tmp += nb_degree_of_freedom;
}
}
++shape_vect;
sh += shapes_size;
}
shape_vect = modified_shapes->begin(nb_degree_of_freedom, lmat_size);
Array<Real>::iterator< Matrix<Real> > lmat = local_mat->begin(lmat_size, lmat_size);
Real * field_val = field_1.values;
for(UInt q = 0; q < vect_size; ++q) {
(*lmat).mul<true, false>(*shape_vect, *shape_vect, *field_val);
++lmat; ++shape_vect; ++field_val;
}
delete modified_shapes;
UInt nb_element = mesh.getNbElement(type, ghost_type);
Array<Real> * int_field_times_shapes = new Array<Real>(nb_element, lmat_size * lmat_size,
- "inte_rho_x_shapes");
+ "inte_rho_x_shapes");
integrator.template integrate<type>(*local_mat, *int_field_times_shapes,
lmat_size * lmat_size, ghost_type, empty_filter);
delete local_mat;
assembleMatrix(*int_field_times_shapes, matrix, nb_degree_of_freedom, type, ghost_type);
delete int_field_times_shapes;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
inline void FEMTemplate<I, S, kind>::inverseMap(const Vector<Real> & real_coords,
UInt element,
const ElementType & type,
Vector<Real> & natural_coords,
const GhostType & ghost_type) const{
AKANTU_DEBUG_IN();
#define INVERSE_MAP(type) \
shape_functions.template inverseMap<type>(real_coords, element, natural_coords, ghost_type);
AKANTU_BOOST_ELEMENT_SWITCH(INVERSE_MAP, AKANTU_NOT_STRUCTURAL_ELEMENT_TYPE);
#undef INVERSE_MAP
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
inline bool FEMTemplate<I, S, kind>::contains(const Vector<Real> & real_coords,
- UInt element,
- const ElementType & type,
- const GhostType & ghost_type) const{
+ UInt element,
+ const ElementType & type,
+ const GhostType & ghost_type) const{
AKANTU_DEBUG_IN();
bool contain = false;
#define CONTAINS(type) \
contain = shape_functions.template contains<type>(real_coords, element, ghost_type);
AKANTU_BOOST_ELEMENT_SWITCH(CONTAINS, AKANTU_NOT_STRUCTURAL_ELEMENT_TYPE);
#undef CONTAINS
AKANTU_DEBUG_OUT();
return contain;
}
/* -------------------------------------------------------------------------- */
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
inline void FEMTemplate<I, S, kind>::computeShapes(const Vector<Real> & real_coords,
- UInt element,
- const ElementType & type,
- Vector<Real> & shapes,
- const GhostType & ghost_type) const{
+ UInt element,
+ const ElementType & type,
+ Vector<Real> & shapes,
+ const GhostType & ghost_type) const{
AKANTU_DEBUG_IN();
#define COMPUTE_SHAPES(type) \
shape_functions.template computeShapes<type>(real_coords,element,shapes,ghost_type);
AKANTU_BOOST_ALL_ELEMENT_SWITCH(COMPUTE_SHAPES);
#undef COMPUTE_SHAPES
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
inline UInt FEMTemplate<I, S, kind>::getNbQuadraturePoints(const ElementType & type,
const GhostType & ghost_type) const {
AKANTU_DEBUG_IN();
UInt nb_quad_points = 0;
#define GET_NB_QUAD(type) \
nb_quad_points = \
integrator. template getQuadraturePoints<type>(ghost_type).cols();
AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_QUAD);
#undef GET_NB_QUAD
AKANTU_DEBUG_OUT();
return nb_quad_points;
}
/* -------------------------------------------------------------------------- */
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
inline const Array<Real> & FEMTemplate<I, S, kind>::getShapes(const ElementType & type,
- const GhostType & ghost_type) const {
+ const GhostType & ghost_type) const {
AKANTU_DEBUG_IN();
const Array<Real> * ret = NULL;
#define GET_SHAPES(type) \
ret = &(shape_functions.getShapes(type, ghost_type));
AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPES);
#undef GET_SHAPES
AKANTU_DEBUG_OUT();
return *ret;
}
/* -------------------------------------------------------------------------- */
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
inline const Array<Real> & FEMTemplate<I, S, kind>::getShapesDerivatives(const ElementType & type,
- const GhostType & ghost_type,
- __attribute__((unused)) UInt id) const {
+ const GhostType & ghost_type,
+ __attribute__((unused)) UInt id) const {
AKANTU_DEBUG_IN();
const Array<Real> * ret = NULL;
#define GET_SHAPES(type) \
ret = &(shape_functions.getShapesDerivatives(type, ghost_type));
AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SHAPES);
#undef GET_SHAPES
AKANTU_DEBUG_OUT();
return *ret;
}
/* -------------------------------------------------------------------------- */
template<template <ElementKind> class I,
template <ElementKind> class S,
ElementKind kind>
inline const Matrix<Real> &
FEMTemplate<I, S, kind>::getQuadraturePoints(const ElementType & type,
const GhostType & ghost_type) const {
AKANTU_DEBUG_IN();
const Matrix<Real> * ret = NULL;
#define GET_QUADS(type) \
ret = &(integrator. template getQuadraturePoints<type>(ghost_type));
AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_QUADS);
#undef GET_QUADS
AKANTU_DEBUG_OUT();
return *ret;
}
/* -------------------------------------------------------------------------- */
__END_AKANTU__
#include "shape_lagrange.hh"
#include "integrator_gauss.hh"
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
template <>
template <>
inline void FEMTemplate<IntegratorGauss, ShapeLagrange, _ek_regular>::
assembleLumpedTemplate<_triangle_6>(const Array<Real> & field_1,
UInt nb_degree_of_freedom,
Array<Real> & lumped,
const Array<Int> & equation_number,
const GhostType & ghost_type) const {
assembleLumpedDiagonalScaling<_triangle_6>(field_1, nb_degree_of_freedom,lumped, equation_number,ghost_type);
}
/* -------------------------------------------------------------------------- */
template <>
template <>
inline void FEMTemplate<IntegratorGauss, ShapeLagrange, _ek_regular>::
assembleLumpedTemplate<_tetrahedron_10>(const Array<Real> & field_1,
UInt nb_degree_of_freedom,
Array<Real> & lumped,
const Array<Int> & equation_number,
const GhostType & ghost_type) const {
assembleLumpedDiagonalScaling<_tetrahedron_10>(field_1, nb_degree_of_freedom,lumped, equation_number,ghost_type);
}
/* -------------------------------------------------------------------------- */
template <>
template <>
inline void
FEMTemplate<IntegratorGauss, ShapeLagrange, _ek_regular>::
assembleLumpedTemplate<_quadrangle_8>(const Array<Real> & field_1,
UInt nb_degree_of_freedom,
Array<Real> & lumped,
const Array<Int> & equation_number,
const GhostType & ghost_type) const {
assembleLumpedDiagonalScaling<_quadrangle_8>(field_1, nb_degree_of_freedom,lumped, equation_number, ghost_type);
}
/* -------------------------------------------------------------------------- */
template<>
template<>
inline void FEMTemplate<IntegratorGauss, ShapeLagrange, _ek_regular>::computeNormalsOnControlPoints<_point_1>(__attribute__((unused))const Array<Real> & field,
Array<Real> & normal,
const GhostType & ghost_type) const {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(mesh.getSpatialDimension() == 1, "Mesh dimension must be 1 to compute normals on points!");
const ElementType type = _point_1;
UInt spatial_dimension = mesh.getSpatialDimension();
//UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
UInt nb_points = getNbQuadraturePoints(type, ghost_type);
UInt nb_element = mesh.getConnectivity(type, ghost_type).getSize();
normal.resize(nb_element * nb_points);
Array<Real>::iterator< Matrix<Real> > normals_on_quad = normal.begin_reinterpret(spatial_dimension,
- nb_points,
- nb_element);
+ nb_points,
+ nb_element);
const Matrix<Real> quads =
integrator.getQuadraturePoints<type>(ghost_type);
Array< std::vector<Element> > segments = mesh.getElementToSubelement(type, ghost_type);
Array<Real> coords = mesh.getNodes();
const Mesh * mesh_segment;
if (mesh.isMeshFacets())
mesh_segment = &(mesh.getMeshParent());
else
mesh_segment = &mesh;
for (UInt elem = 0; elem < nb_element; ++elem) {
UInt nb_segment = segments(elem).size();
AKANTU_DEBUG_ASSERT(nb_segment > 0, "Impossible to compute a normal on a point connected to 0 segments");
Real normal_value = 1;
if (nb_segment == 1) {
Element segment = segments(elem)[0];
const Array<UInt> & segment_connectivity = mesh_segment->getConnectivity(segment.type, segment.ghost_type);
//const Vector<UInt> & segment_points = segment_connectivity.begin(Mesh::getNbNodesPerElement(segment.type))[segment.element];
Real difference;
if(segment_connectivity(0) == elem) {
difference = coords(elem)-coords(segment_connectivity(1));
} else {
difference = coords(elem)-coords(segment_connectivity(0));
}
normal_value = difference/std::abs(difference);
}
for(UInt n(0); n < nb_points; ++n) {
(*normals_on_quad)(0, n) = normal_value;
}
++normals_on_quad;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/* Shape Linked specialization */
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
__END_AKANTU__
#include "shape_linked.hh"
__BEGIN_AKANTU__
#if defined(AKANTU_STRUCTURAL_MECHANICS)
/* -------------------------------------------------------------------------- */
template <>
inline const Array<Real> &
FEMTemplate<IntegratorGauss, ShapeLinked, _ek_structural>::getShapesDerivatives(const ElementType & type,
const GhostType & ghost_type,
UInt id) const {
AKANTU_DEBUG_IN();
const Array<Real> * ret = NULL;
#define GET_SHAPES(type) \
ret = &(shape_functions.getShapesDerivatives(type, ghost_type, id));
AKANTU_BOOST_STRUCTURAL_ELEMENT_SWITCH(GET_SHAPES);
#undef GET_SHAPES
AKANTU_DEBUG_OUT();
return *ret;
}
#endif
diff --git a/src/fem/group_manager.cc b/src/fem/group_manager.cc
index 96dbfe5f5..645a9dde2 100644
--- a/src/fem/group_manager.cc
+++ b/src/fem/group_manager.cc
@@ -1,800 +1,740 @@
/**
* @file group_manager.cc
*
* @author Dana Christen <dana.christen@gmail.com>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author David Kammer <david.kammer@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date Wed Mar 06 09:30:00 2013
*
* @brief Stores information about ElementGroup and NodeGroup
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "group_manager.hh"
#include "mesh.hh"
#include "aka_csr.hh"
#include "mesh_utils.hh"
#include "element_group.hh"
#include "node_group.hh"
#include "data_accessor.hh"
#include "distributed_synchronizer.hh"
/* -------------------------------------------------------------------------- */
#include <sstream>
#include <algorithm>
#include <iterator>
#include <list>
#include <queue>
#include <numeric>
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
GroupManager::GroupManager(const Mesh & mesh,
const ID & id,
const MemoryID & mem_id) : id(id),
memory_id(mem_id),
mesh(mesh) {
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
GroupManager::~GroupManager() {
ElementGroups::iterator eit = element_groups.begin();
ElementGroups::iterator eend = element_groups.end();
for(; eit != eend ; ++eit) delete (eit->second);
NodeGroups::iterator nit = node_groups.begin();
NodeGroups::iterator nend = node_groups.end();
for(; nit != nend ; ++nit) delete (nit->second);
}
/* -------------------------------------------------------------------------- */
NodeGroup & GroupManager::createNodeGroup(const std::string & group_name,
bool replace_group) {
AKANTU_DEBUG_IN();
NodeGroups::iterator it = node_groups.find(group_name);
if(it != node_groups.end()) {
if (replace_group) {
it->second->empty();
AKANTU_DEBUG_OUT();
return *(it->second);
}
else
AKANTU_EXCEPTION("Trying to create a node group that already exists:" << group_name << "_nodes");
}
NodeGroup * node_group = new NodeGroup(group_name, id + ":" + group_name + "_node_group",
memory_id);
node_groups[group_name] = node_group;
AKANTU_DEBUG_OUT();
return *node_group;
}
/* -------------------------------------------------------------------------- */
void GroupManager::destroyNodeGroup(const std::string & group_name) {
AKANTU_DEBUG_IN();
NodeGroups::iterator nit = node_groups.find(group_name);
NodeGroups::iterator nend = node_groups.end();
if (nit != nend) {
delete (nit->second);
node_groups.erase(nit);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
ElementGroup & GroupManager::createElementGroup(const std::string & group_name,
UInt dimension,
bool replace_group) {
AKANTU_DEBUG_IN();
NodeGroup & new_node_group = createNodeGroup(group_name + "_nodes", replace_group);
ElementGroups::iterator it = element_groups.find(group_name);
if(it != element_groups.end()) {
if (replace_group) {
it->second->empty();
AKANTU_DEBUG_OUT();
return *(it->second);
}
else
AKANTU_EXCEPTION("Trying to create a element group that already exists:" << group_name);
}
ElementGroup * element_group = new ElementGroup(group_name, mesh, new_node_group,
dimension,
id + ":" + group_name + "_element_group",
memory_id);
node_groups[group_name + "_nodes"] = &new_node_group;
element_groups[group_name] = element_group;
AKANTU_DEBUG_OUT();
return *element_group;
}
/* -------------------------------------------------------------------------- */
void GroupManager::destroyElementGroup(const std::string & group_name,
bool destroy_node_group) {
AKANTU_DEBUG_IN();
ElementGroups::iterator eit = element_groups.find(group_name);
ElementGroups::iterator eend = element_groups.end();
if (eit != eend) {
if (destroy_node_group)
destroyNodeGroup(eit->second->getNodeGroup().getName());
delete (eit->second);
element_groups.erase(eit);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
ElementGroup & GroupManager::createElementGroup(const std::string & group_name,
UInt dimension,
NodeGroup & node_group) {
AKANTU_DEBUG_IN();
if(element_groups.find(group_name) != element_groups.end()) {
AKANTU_EXCEPTION("Trying to create a element group that already exists:" << group_name);
}
ElementGroup * element_group = new ElementGroup(group_name, mesh, node_group,
dimension,
id + ":" + group_name + "_element_group",
memory_id);
element_groups[group_name] = element_group;
AKANTU_DEBUG_OUT();
return *element_group;
}
/* -------------------------------------------------------------------------- */
class ClusterSynchronizer : public DataAccessor {
typedef std::set< std::pair<UInt, UInt> > DistantIDs;
public:
ClusterSynchronizer(GroupManager & group_manager,
UInt element_dimension,
std::string cluster_name_prefix,
ByElementTypeUInt & element_to_fragment,
DistributedSynchronizer & distributed_synchronizer,
UInt nb_cluster) :
group_manager(group_manager),
element_dimension(element_dimension),
cluster_name_prefix(cluster_name_prefix),
element_to_fragment(element_to_fragment),
distributed_synchronizer(distributed_synchronizer),
nb_cluster(nb_cluster) { }
UInt synchronize() {
StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator();
UInt rank = comm.whoAmI();
UInt nb_proc = comm.getNbProc();
/// find starting index to renumber local clusters
Array<UInt> nb_cluster_per_proc(nb_proc);
nb_cluster_per_proc(rank) = nb_cluster;
comm.allGather(nb_cluster_per_proc.storage(), 1);
starting_index = std::accumulate(nb_cluster_per_proc.begin(),
nb_cluster_per_proc.begin() + rank,
0);
UInt global_nb_fragment = std::accumulate(nb_cluster_per_proc.begin() + rank,
nb_cluster_per_proc.end(),
starting_index);
/// create the local to distant cluster pairs with neighbors
distributed_synchronizer.computeBufferSize(*this, _gst_gm_clusters);
distributed_synchronizer.asynchronousSynchronize(*this, _gst_gm_clusters);
distributed_synchronizer.waitEndSynchronize(*this, _gst_gm_clusters);
/// count total number of pairs
Array<Int> nb_pairs(nb_proc);
nb_pairs(rank) = distant_ids.size();
comm.allGather(nb_pairs.storage(), 1);
UInt total_nb_pairs = std::accumulate(nb_pairs.begin(), nb_pairs.end(), 0);
/// generate pairs global array
UInt local_pair_index = std::accumulate(nb_pairs.storage(),
nb_pairs.storage() + rank, 0);
Array<UInt> total_pairs(total_nb_pairs, 2);
DistantIDs::iterator ids_it = distant_ids.begin();
DistantIDs::iterator ids_end = distant_ids.end();
for (; ids_it != ids_end; ++ids_it, ++local_pair_index) {
total_pairs(local_pair_index, 0) = ids_it->first;
total_pairs(local_pair_index, 1) = ids_it->second;
}
/// communicate pairs to all processors
nb_pairs *= 2;
comm.allGatherV(total_pairs.storage(), nb_pairs.storage());
/// renumber clusters
Array<UInt>::iterator<Vector<UInt> > pairs_it = total_pairs.begin(2);
Array<UInt>::iterator<Vector<UInt> > pairs_end = total_pairs.end(2);
/// generate fragment list
std::vector< std::set<UInt> > global_clusters;
UInt total_nb_cluster = 0;
Array<bool> is_fragment_in_cluster(global_nb_fragment, 1, false);
std::queue<UInt> fragment_check_list;
while (total_pairs.getSize() != 0) {
/// create a new cluster
++total_nb_cluster;
global_clusters.resize(total_nb_cluster);
std::set<UInt> & current_cluster = global_clusters[total_nb_cluster - 1];
UInt first_fragment = total_pairs(0, 0);
UInt second_fragment = total_pairs(0, 1);
total_pairs.erase(0);
fragment_check_list.push(first_fragment);
fragment_check_list.push(second_fragment);
while (!fragment_check_list.empty()) {
UInt current_fragment = fragment_check_list.front();
UInt * total_pairs_end = total_pairs.storage() + total_pairs.getSize() * 2;
UInt * fragment_found = std::find(total_pairs.storage(),
total_pairs_end,
current_fragment);
if (fragment_found != total_pairs_end) {
UInt position = fragment_found - total_pairs.storage();
UInt pair = position / 2;
UInt other_index = (position + 1) % 2;
fragment_check_list.push(total_pairs(pair, other_index));
total_pairs.erase(pair);
}
else {
fragment_check_list.pop();
current_cluster.insert(current_fragment);
is_fragment_in_cluster(current_fragment) = true;
}
}
}
/// add to FragmentToCluster all local fragments
for (UInt c = 0; c < global_nb_fragment; ++c) {
if (!is_fragment_in_cluster(c)) {
++total_nb_cluster;
global_clusters.resize(total_nb_cluster);
std::set<UInt> & current_cluster = global_clusters[total_nb_cluster - 1];
current_cluster.insert(c);
}
}
/// reorganize element groups to match global clusters
for (UInt c = 0; c < global_clusters.size(); ++c) {
/// create new element group corresponding to current cluster
std::stringstream sstr;
sstr << cluster_name_prefix << "_" << c;
ElementGroup & cluster = group_manager.createElementGroup(sstr.str(),
element_dimension,
true);
std::set<UInt>::iterator it = global_clusters[c].begin();
std::set<UInt>::iterator end = global_clusters[c].end();
/// append to current element group all fragments that belong to
/// the same cluster if they exist
for (; it != end; ++it) {
Int local_index = *it - starting_index;
if (local_index < 0 || local_index >= Int(nb_cluster)) continue;
std::stringstream tmp_sstr;
tmp_sstr << "tmp_" << cluster_name_prefix << "_" << local_index;
GroupManager::element_group_iterator eg_it
= group_manager.element_group_find(tmp_sstr.str());
AKANTU_DEBUG_ASSERT(eg_it != group_manager.element_group_end(),
"Temporary fragment \""<< tmp_sstr.str() << "\" not found");
cluster.append(*(eg_it->second));
group_manager.destroyElementGroup(tmp_sstr.str(), true);
}
}
return total_nb_cluster;
}
private:
/// functions for parallel communications
inline UInt getNbDataForElements(const Array<Element> & elements,
SynchronizationTag tag) const {
if (tag == _gst_gm_clusters)
return elements.getSize() * sizeof(UInt);
return 0;
}
inline void packElementData(CommunicationBuffer & buffer,
const Array<Element> & elements,
SynchronizationTag tag) const {
if (tag != _gst_gm_clusters) return;
Array<Element>::const_iterator<> el_it = elements.begin();
Array<Element>::const_iterator<> el_end = elements.end();
for (; el_it != el_end; ++el_it) {
const Element & el = *el_it;
/// for each element pack its global cluster index
buffer << element_to_fragment(el.type, el.ghost_type)(el.element) + starting_index;
}
}
inline void unpackElementData(CommunicationBuffer & buffer,
const Array<Element> & elements,
SynchronizationTag tag) {
if (tag != _gst_gm_clusters) return;
Array<Element>::const_iterator<> el_it = elements.begin();
Array<Element>::const_iterator<> el_end = elements.end();
for (; el_it != el_end; ++el_it) {
UInt distant_cluster;
buffer >> distant_cluster;
const Element & el = *el_it;
UInt local_cluster = element_to_fragment(el.type, el.ghost_type)(el.element) + starting_index;
distant_ids.insert(std::make_pair(local_cluster, distant_cluster));
}
}
private:
GroupManager & group_manager;
UInt element_dimension;
std::string cluster_name_prefix;
ByElementTypeUInt & element_to_fragment;
DistributedSynchronizer & distributed_synchronizer;
UInt nb_cluster;
DistantIDs distant_ids;
UInt starting_index;
};
/* -------------------------------------------------------------------------- */
/// \todo this function doesn't work in 1D
UInt GroupManager::createBoundaryGroupFromGeometry() {
UInt spatial_dimension = mesh.getSpatialDimension();
return createClusters(spatial_dimension - 1, "boundary");
}
/* -------------------------------------------------------------------------- */
//// \todo if needed element list construction can be optimized by
//// templating the filter class
UInt GroupManager::createClusters(UInt element_dimension,
std::string cluster_name_prefix,
const GroupManager::ClusteringFilter & filter,
DistributedSynchronizer * distributed_synchronizer,
Mesh * mesh_facets) {
AKANTU_DEBUG_IN();
UInt nb_proc = StaticCommunicator::getStaticCommunicator().getNbProc();
std::string tmp_cluster_name_prefix = cluster_name_prefix;
ByElementTypeUInt * element_to_fragment = NULL;
if (nb_proc > 1 && distributed_synchronizer) {
element_to_fragment = new ByElementTypeUInt;
mesh.initByElementTypeArray(*element_to_fragment, 1, element_dimension,
false, _ek_not_defined, true);
tmp_cluster_name_prefix = "tmp_" + tmp_cluster_name_prefix;
}
/// Get facets
bool mesh_facets_created = false;
if (!mesh_facets && element_dimension > 0) {
mesh_facets = new Mesh(mesh.getSpatialDimension(),
mesh.getNodes().getID(),
"mesh_facets_for_clusters");
mesh_facets->defineMeshParent(mesh);
MeshUtils::buildAllFacets(mesh, *mesh_facets,
element_dimension - 1,
distributed_synchronizer);
}
std::list<Element> unseen_elements;
Element el;
for (ghost_type_t::iterator gt = ghost_type_t::begin();
gt != ghost_type_t::end(); ++gt) {
GhostType ghost_type = *gt;
el.ghost_type = ghost_type;
Mesh::type_iterator type_it = mesh.firstType(element_dimension,
ghost_type, _ek_not_defined);
Mesh::type_iterator type_end = mesh.lastType (element_dimension,
ghost_type, _ek_not_defined);
for (; type_it != type_end; ++type_it) {
el.type = *type_it;
el.kind = Mesh::getKind(*type_it);
UInt nb_element = mesh.getNbElement(*type_it, ghost_type);
for (UInt e = 0; e < nb_element; ++e) {
el.element = e;
if (filter(el))
unseen_elements.push_back(el);
}
}
}
Array<bool> checked_node(mesh.getNbNodes(), 1, false);
UInt nb_cluster = 0;
/// keep looping until all elements are seen
while(!unseen_elements.empty()) {
/// create a new cluster
std::stringstream sstr;
sstr << tmp_cluster_name_prefix << "_" << nb_cluster;
ElementGroup & cluster = createElementGroup(sstr.str(),
element_dimension,
true);
++nb_cluster;
/// initialize the queue of elements to check in the current cluster
std::list<Element>::iterator uns_it = unseen_elements.begin();
Element uns_el = *uns_it;
unseen_elements.erase(uns_it);
// point element are cluster by themself
if(element_dimension == 0) {
cluster.add(uns_el);
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(el.type);
Vector<UInt> connect =
mesh.getConnectivity(uns_el.type, uns_el.ghost_type).begin(nb_nodes_per_element)[uns_el.element];
for (UInt n = 0; n < nb_nodes_per_element; ++n) {
/// add element's nodes to the cluster
UInt node = connect[n];
if (!checked_node(node)) {
cluster.addNode(node);
checked_node(node) = true;
}
}
continue;
}
std::queue<Element> element_to_add;
element_to_add.push(uns_el);
/// keep looping until current cluster is complete (no more
/// connected elements)
while(!element_to_add.empty()) {
/// take first element and erase it in the queue
Element el = element_to_add.front();
element_to_add.pop();
/// if parallel, store cluster index per element
if (nb_proc > 1 && distributed_synchronizer)
(*element_to_fragment)(el.type, el.ghost_type)(el.element) = nb_cluster - 1;
/// add current element to the cluster
cluster.add(el);
const Array<Element> & element_to_facet
= mesh_facets->getSubelementToElement(el.type, el.ghost_type);
UInt nb_facet_per_element = element_to_facet.getNbComponent();
for (UInt f = 0; f < nb_facet_per_element; ++f) {
const Element & facet = element_to_facet(el.element, f);
if (facet == ElementNull) continue;
const std::vector<Element> & connected_elements
= mesh_facets->getElementToSubelement(facet.type, facet.ghost_type)(facet.element);
for (UInt elem = 0; elem < connected_elements.size(); ++elem) {
const Element & check_el = connected_elements[elem];
if (check_el == ElementNull || check_el == el) continue;
std::list<Element>::iterator it_clus = std::find(unseen_elements.begin(),
unseen_elements.end(),
check_el);
/// if neighbor not seen yet, add it to check list and
/// remove it from unseen elements
if(it_clus != unseen_elements.end()) {
unseen_elements.erase(it_clus);
element_to_add.push(check_el);
}
}
}
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(el.type);
Vector<UInt> connect =
mesh.getConnectivity(el.type, el.ghost_type).begin(nb_nodes_per_element)[el.element];
for (UInt n = 0; n < nb_nodes_per_element; ++n) {
/// add element's nodes to the cluster
UInt node = connect[n];
if (!checked_node(node)) {
cluster.addNode(node);
checked_node(node) = true;
}
}
}
}
if (nb_proc > 1 && distributed_synchronizer) {
ClusterSynchronizer cluster_synchronizer(*this, element_dimension,
cluster_name_prefix,
*element_to_fragment,
*distributed_synchronizer,
nb_cluster);
nb_cluster = cluster_synchronizer.synchronize();
delete element_to_fragment;
}
if (mesh_facets_created)
delete mesh_facets;
AKANTU_DEBUG_OUT();
return nb_cluster;
}
-/* -------------------------------------------------------------------------- */
-void GroupManager::createMeshDataFromClusters(const std::string cluster_name_prefix) {
- AKANTU_DEBUG_IN();
-
- std::string mesh_data_name(cluster_name_prefix);
- Mesh & mesh_not_const = const_cast<Mesh &>(mesh);
-
- /// loop over clusters
- for(const_element_group_iterator it(element_group_begin());
- it != element_group_end(); ++it) {
-
- /// search for prefix in cluster name
- size_t pos = it->first.find(cluster_name_prefix);
-
- /// if not found continue
- if (pos == std::string::npos) continue;
-
- AKANTU_DEBUG_ASSERT(pos == 0, "cluster_name_prefix is wrong");
-
- /// get cluster index
- std::string cluster_index_string
- = it->first.substr(cluster_name_prefix.size() + 1);
- std::stringstream sstr(cluster_index_string.c_str());
- UInt cluster_index;
- sstr >> cluster_index;
-
- AKANTU_DEBUG_ASSERT(!sstr.fail(), "cluster_index is not an integer");
-
- /// loop over cluster types
- for (ghost_type_t::iterator gt = ghost_type_t::begin();
- gt != ghost_type_t::end(); ++gt) {
- GhostType ghost_type = *gt;
-
- ElementGroup::type_iterator type_it
- = it->second->firstType(_all_dimensions, ghost_type);
- ElementGroup::type_iterator type_end
- = it->second->lastType(_all_dimensions, ghost_type);
-
- for(; type_it != type_end; ++type_it) {
- ElementType type = *type_it;
-
- /// init mesh data
- Array<UInt> * mesh_data
- = mesh_not_const.getDataPointer<UInt>(mesh_data_name, type, ghost_type, 1);
-
- ElementGroup::const_element_iterator el_it
- = it->second->element_begin(type, ghost_type);
- ElementGroup::const_element_iterator el_end
- = it->second->element_end(type, ghost_type);
-
- /// fill mesh data with cluster index
- for (; el_it != el_end; ++el_it)
- (*mesh_data)(*el_it) = cluster_index;
- }
- }
- }
-
- AKANTU_DEBUG_OUT();
-}
-
/* -------------------------------------------------------------------------- */
template<typename T>
void GroupManager::createGroupsFromMeshData(const std::string & dataset_name) {
std::set<std::string> group_names;
const ByElementTypeArray<T> & datas = mesh.getData<T>(dataset_name);
typedef typename ByElementTypeArray<T>::type_iterator type_iterator;
for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) {
type_iterator type_it = datas.firstType(_all_dimensions, *gt);
type_iterator type_end = datas.lastType(_all_dimensions, *gt);
for (; type_it != type_end; ++type_it) {
const Array<T> & dataset = datas(*type_it, *gt);
UInt nb_element = mesh.getNbElement(*type_it, *gt);
AKANTU_DEBUG_ASSERT(dataset.getSize() == nb_element,
"Not the same number of elements in the map from MeshData and in the mesh!");
for(UInt e(0); e < nb_element; ++e) {
std::stringstream sstr; sstr << dataset(e);
group_names.insert(sstr.str());
}
}
}
/// @todo sharing the group names in parallel to avoid trouble with groups
/// present only on a sub set of processors
std::set<std::string>::iterator git = group_names.begin();
std::set<std::string>::iterator gend = group_names.end();
for (;git != gend; ++git) createElementGroup(*git, _all_dimensions);
Element el;
for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) {
el.ghost_type = *gt;
type_iterator type_it = datas.firstType(_all_dimensions, *gt);
type_iterator type_end = datas.lastType(_all_dimensions, *gt);
for (; type_it != type_end; ++type_it) {
el.type = *type_it;
const Array<T> & dataset = datas(*type_it, *gt);
UInt nb_element = mesh.getNbElement(*type_it, *gt);
AKANTU_DEBUG_ASSERT(dataset.getSize() == nb_element,
"Not the same number of elements in the map from MeshData and in the mesh!");
UInt nb_nodes_per_element = mesh.getNbNodesPerElement(el.type);
Array<UInt>::const_iterator< Vector<UInt> > cit =
mesh.getConnectivity(*type_it, *gt).begin(nb_nodes_per_element);
for(UInt e(0); e < nb_element; ++e, ++cit) {
el.element = e;
std::stringstream sstr; sstr << dataset(e);
ElementGroup & group = getElementGroup(sstr.str());
group.add(el);
const Vector<UInt> & connect = *cit;
for (UInt n = 0; n < nb_nodes_per_element; ++n) {
UInt node = connect[n];
group.addNode(node);
}
}
}
}
git = group_names.begin();
for (;git != gend; ++git) getElementGroup(*git).getNodeGroup().removeDuplicate();
}
template void GroupManager::createGroupsFromMeshData<std::string>(const std::string & dataset_name);
template void GroupManager::createGroupsFromMeshData<UInt>(const std::string & dataset_name);
/* -------------------------------------------------------------------------- */
void GroupManager::createElementGroupFromNodeGroup(const std::string & name,
const std::string & node_group_name,
UInt dimension) {
NodeGroup & node_group = getNodeGroup(node_group_name);
ElementGroup & group = createElementGroup(name, dimension, node_group);
CSR<Element> node_to_elem;
MeshUtils::buildNode2Elements(mesh, node_to_elem, _all_dimensions);
std::set<Element> seen;
Array<UInt>::const_iterator<> itn = node_group.begin();
Array<UInt>::const_iterator<> endn = node_group.end();
for (;itn != endn; ++itn) {
CSR<Element>::iterator ite = node_to_elem.begin(*itn);
CSR<Element>::iterator ende = node_to_elem.end(*itn);
for (;ite != ende; ++ite) {
const Element & elem = *ite;
if(dimension != _all_dimensions && dimension != Mesh::getSpatialDimension(elem.type)) continue;
if(seen.find(elem) != seen.end()) continue;
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(elem.type);
Array<UInt>::const_iterator< Vector<UInt> > conn_it =
mesh.getConnectivity(elem.type, elem.ghost_type).begin(nb_nodes_per_element);
const Vector<UInt> & conn = conn_it[elem.element];
UInt count = 0;
for (UInt n = 0; n < conn.size(); ++n) {
count += (node_group.getNodes().find(conn(n)) != -1 ? 1 : 0);
}
if(count == nb_nodes_per_element) group.add(elem);
seen.insert(elem);
}
}
}
/* -------------------------------------------------------------------------- */
void GroupManager::printself(std::ostream & stream, int indent) const {
std::string space;
for(Int i = 0; i < indent; i++, space += AKANTU_INDENT);
stream << space << "GroupManager [" << std::endl;
std::set<std::string> node_group_seen;
for(const_element_group_iterator it(element_group_begin()); it != element_group_end(); ++it) {
it->second->printself(stream, indent + 1);
node_group_seen.insert(it->second->getNodeGroup().getName());
}
for(const_node_group_iterator it(node_group_begin()); it != node_group_end(); ++it) {
if(node_group_seen.find(it->second->getName()) == node_group_seen.end())
it->second->printself(stream, indent + 1);
}
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
UInt GroupManager::getNbElementGroups(UInt dimension) const {
if(dimension == _all_dimensions) return element_groups.size();
ElementGroups::const_iterator it = element_groups.begin();
ElementGroups::const_iterator end = element_groups.end();
UInt count = 0;
for(;it != end; ++it) count += (it->second->getDimension() == dimension);
return count;
}
__END_AKANTU__
diff --git a/src/fem/group_manager.hh b/src/fem/group_manager.hh
index ac3e2ab97..db2409983 100644
--- a/src/fem/group_manager.hh
+++ b/src/fem/group_manager.hh
@@ -1,207 +1,204 @@
/**
* @file group_manager.hh
*
* @author Dana Christen <dana.christen@gmail.com>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author David Kammer <david.kammer@epfl.ch>
*
* @date Wed Mar 06 09:30:00 2013
*
* @brief Stores information relevent to the notion of element and nodes groups.
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_GROUP_MANAGER_HH__
#define __AKANTU_GROUP_MANAGER_HH__
#include <set>
#include "aka_common.hh"
#include "by_element_type.hh"
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
class ElementGroup;
class NodeGroup;
class Mesh;
class Element;
class DistributedSynchronizer;
class GroupManager {
/* ------------------------------------------------------------------------ */
/* Typedefs */
/* ------------------------------------------------------------------------ */
private:
typedef std::map<std::string, ElementGroup *> ElementGroups;
typedef std::map<std::string, NodeGroup *> NodeGroups;
public:
typedef std::set<ElementType> GroupManagerTypeSet;
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
GroupManager(const Mesh & mesh,
const ID & id = "group_manager",
const MemoryID & memory_id = 0);
- ~GroupManager();
+ virtual ~GroupManager();
/* ------------------------------------------------------------------------ */
/* Groups iterators */
/* ------------------------------------------------------------------------ */
- public:
+public:
typedef NodeGroups::iterator node_group_iterator;
typedef ElementGroups::iterator element_group_iterator;
typedef NodeGroups::const_iterator const_node_group_iterator;
typedef ElementGroups::const_iterator const_element_group_iterator;
#define AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION(group_type, \
function, \
param_in, \
param_out) \
inline BOOST_PP_CAT(BOOST_PP_CAT(const_, group_type), _iterator) \
BOOST_PP_CAT(BOOST_PP_CAT(group_type, _), function)(param_in) const { \
return BOOST_PP_CAT(group_type, s).function(param_out); \
}; \
\
inline BOOST_PP_CAT(group_type, _iterator) \
BOOST_PP_CAT(BOOST_PP_CAT(group_type, _), function)(param_in) { \
return BOOST_PP_CAT(group_type, s).function(param_out); \
}
#define AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(group_type, function) \
AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION(group_type, function, BOOST_PP_EMPTY(), BOOST_PP_EMPTY())
AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(node_group, begin);
AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(node_group, end );
AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(element_group, begin);
AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION_NP(element_group, end );
AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION(element_group, find, const std::string & name, name);
AKANTU_GROUP_MANAGER_DEFINE_ITERATOR_FUNCTION(node_group, find, const std::string & name, name);
/* ------------------------------------------------------------------------ */
/* Clustering filter */
/* ------------------------------------------------------------------------ */
public:
class ClusteringFilter {
public:
virtual bool operator() (const Element &) const {
return true;
}
};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// create a node group
NodeGroup & createNodeGroup(const std::string & group_name,
bool replace_group = false);
/// destroy a node group
void destroyNodeGroup(const std::string & group_name);
/// create an element group and the associated node group
ElementGroup & createElementGroup(const std::string & group_name,
UInt dimension,
bool replace_group = false);
/// destroy an element group and the associated node group
void destroyElementGroup(const std::string & group_name,
bool destroy_node_group = false);
/// create a element group using an existing node group
ElementGroup & createElementGroup(const std::string & group_name, UInt dimension, NodeGroup & node_group);
- /// create mesh data based on clusters
- void createMeshDataFromClusters(const std::string cluster_name_prefix);
-
/// create groups based on values stored in a given mesh data
template <typename T>
void createGroupsFromMeshData(const std::string & dataset_name);
/// create boundaries group by a clustering algorithm \todo extend to parallel
UInt createBoundaryGroupFromGeometry();
/// create element clusters for a given dimension
UInt createClusters(UInt element_dimension,
- std::string cluster_name_prefix = "cluster_",
+ std::string cluster_name_prefix = "cluster",
const ClusteringFilter & filter = ClusteringFilter(),
DistributedSynchronizer * distributed_synchronizer = NULL,
Mesh * mesh_facets = NULL);
/// Create an ElementGroup based on a NodeGroup
void createElementGroupFromNodeGroup(const std::string & name,
const std::string & node_group,
UInt dimension = _all_dimensions);
- void printself(std::ostream & stream, int indent = 0) const;
+ virtual void printself(std::ostream & stream, int indent = 0) const;
/* ------------------------------------------------------------------------ */
/* Accessor */
/* ------------------------------------------------------------------------ */
public:
inline const ElementGroup & getElementGroup(const std::string & name) const;
inline const NodeGroup & getNodeGroup(const std::string & name) const;
inline ElementGroup & getElementGroup(const std::string & name);
inline NodeGroup & getNodeGroup(const std::string & name);
UInt getNbElementGroups(UInt dimension = _all_dimensions) const;
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
-private:
+protected:
/// id to create element and node groups
ID id;
/// memory_id to create element and node groups
MemoryID memory_id;
/// list of the node groups managed
NodeGroups node_groups;
/// list of the node groups managed
ElementGroups element_groups;
/// Mesh to which the element belongs
const Mesh & mesh;
};
/// standard output stream operator
inline std::ostream & operator <<(std::ostream & stream, const GroupManager & _this)
{
_this.printself(stream);
return stream;
}
__END_AKANTU__
#include "group_manager_inline_impl.cc"
#endif /* __AKANTU_GROUP_MANAGER_HH__ */
diff --git a/src/model/solid_mechanics/fragment_manager.cc b/src/model/solid_mechanics/fragment_manager.cc
new file mode 100644
index 000000000..f62d5540e
--- /dev/null
+++ b/src/model/solid_mechanics/fragment_manager.cc
@@ -0,0 +1,577 @@
+/**
+ * @file fragment_manager.cc
+ * @author Marco Vocialta <marco.vocialta@epfl.ch>
+ * @date Mon Jan 20 14:33:49 2014
+ *
+ * @brief Group manager to handle fragments
+ *
+ * @section LICENSE
+ *
+ * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
+ * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
+ *
+ * Akantu is free software: you can redistribute it and/or modify it under the
+ * terms of the GNU Lesser General Public License as published by the Free
+ * Software Foundation, either version 3 of the License, or (at your option) any
+ * later version.
+ *
+ * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
+ * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+ * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+ * details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+/* -------------------------------------------------------------------------- */
+#include "fragment_manager.hh"
+#include "material_cohesive.hh"
+#include <numeric>
+
+__BEGIN_AKANTU__
+
+/* -------------------------------------------------------------------------- */
+FragmentManager::FragmentManager(SolidMechanicsModelCohesive & model,
+ bool dump_data,
+ const ID & id,
+ const MemoryID & memory_id) :
+ GroupManager(model.getMesh(), id, memory_id),
+ model(model),
+ mass_center(0, model.getSpatialDimension(), "mass_center"),
+ mass(0, model.getSpatialDimension(), "mass"),
+ velocity(0, model.getSpatialDimension(), "velocity"),
+ inertia_moments(0, model.getSpatialDimension(), "inertia_moments"),
+ quad_coordinates("quad_coordinates", id),
+ mass_density("mass_density", id),
+ dump_data(dump_data) {
+ AKANTU_DEBUG_IN();
+
+ UInt spatial_dimension = mesh.getSpatialDimension();
+
+ /// compute quadrature points' coordinates
+ mesh.initByElementTypeArray(quad_coordinates,
+ spatial_dimension,
+ spatial_dimension,
+ _not_ghost);
+
+ model.getFEM().interpolateOnQuadraturePoints(model.getMesh().getNodes(),
+ quad_coordinates);
+
+ /// store mass density per quadrature point
+ mesh.initByElementTypeArray(mass_density,
+ 1,
+ spatial_dimension,
+ _not_ghost);
+
+ storeMassDensityPerQuadraturePoint();
+
+ AKANTU_DEBUG_OUT();
+}
+
+/* -------------------------------------------------------------------------- */
+class CohesiveElementFilter : public GroupManager::ClusteringFilter {
+public:
+ CohesiveElementFilter(const SolidMechanicsModelCohesive & model,
+ const Real max_damage = 1.) :
+ model(model), is_unbroken(max_damage) {}
+
+ bool operator() (const Element & el) const {
+ if (el.kind == _ek_regular)
+ return true;
+
+ const Array<UInt> & el_id_by_mat = model.getElementIndexByMaterial(el.type,
+ el.ghost_type);
+
+ const MaterialCohesive & mat
+ = static_cast<const MaterialCohesive &>
+ (model.getMaterial(el_id_by_mat(el.element, 0)));
+
+ UInt el_index = el_id_by_mat(el.element, 1);
+ UInt nb_quad_per_element
+ = model.getFEM("CohesiveFEM").getNbQuadraturePoints(el.type, el.ghost_type);
+
+ Vector<Real> damage
+ = mat.getDamage(el.type, el.ghost_type).begin(nb_quad_per_element)[el_index];
+
+ UInt unbroken_quads = std::count_if(damage.storage(),
+ damage.storage() + nb_quad_per_element,
+ is_unbroken);
+
+ if (unbroken_quads > 0)
+ return true;
+ return false;
+ }
+
+private:
+
+ struct IsUnbrokenFunctor {
+ IsUnbrokenFunctor(const Real & max_damage) : max_damage(max_damage) {}
+ bool operator() (const Real & x) {return x < max_damage;}
+ const Real max_damage;
+ };
+
+ const SolidMechanicsModelCohesive & model;
+ const IsUnbrokenFunctor is_unbroken;
+};
+
+/* -------------------------------------------------------------------------- */
+void FragmentManager::buildFragments() {
+ AKANTU_DEBUG_IN();
+
+#if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT)
+ DistributedSynchronizer * cohesive_synchronizer
+ = const_cast<DistributedSynchronizer *>(model.getCohesiveSynchronizer());
+
+ if (cohesive_synchronizer) {
+ cohesive_synchronizer->computeBufferSize(model, _gst_smmc_damage);
+ cohesive_synchronizer->asynchronousSynchronize(model, _gst_smmc_damage);
+ cohesive_synchronizer->waitEndSynchronize(model, _gst_smmc_damage);
+ }
+#endif
+
+ DistributedSynchronizer & synchronizer
+ = const_cast<DistributedSynchronizer &>(model.getSynchronizer());
+
+ Mesh & mesh_facets = const_cast<Mesh &>(model.getMeshFacets());
+
+ UInt spatial_dimension = model.getSpatialDimension();
+ std::string fragment_prefix("fragment");
+
+ /// generate fragments
+ global_nb_fragment = createClusters(spatial_dimension,
+ fragment_prefix,
+ CohesiveElementFilter(model),
+ &synchronizer,
+ &mesh_facets);
+
+ nb_fragment = getNbElementGroups(spatial_dimension);
+ fragment_index.resize(nb_fragment);
+
+ UInt * fragment_index_it = fragment_index.storage();
+
+ /// loop over fragments
+ for(const_element_group_iterator it(element_group_begin());
+ it != element_group_end(); ++it, ++fragment_index_it) {
+
+ /// get fragment index
+ std::string fragment_index_string
+ = it->first.substr(fragment_prefix.size() + 1);
+ std::stringstream sstr(fragment_index_string.c_str());
+ sstr >> *fragment_index_it;
+
+ AKANTU_DEBUG_ASSERT(!sstr.fail(), "fragment_index is not an integer");
+ }
+
+ /// compute fragments' mass
+ computeMass();
+
+ if (dump_data) {
+ createDumpDataArray(fragment_index, "fragments");
+ createDumpDataArray(mass, "fragments mass");
+ }
+
+ AKANTU_DEBUG_OUT();
+}
+
+/* -------------------------------------------------------------------------- */
+void FragmentManager::computeMass() {
+ AKANTU_DEBUG_IN();
+
+ UInt spatial_dimension = model.getSpatialDimension();
+
+ /// create a unit field per quadrature point, since to compute mass
+ /// it's neccessary to integrate only density
+ ByElementTypeReal unit_field("unit_field", id);
+ mesh.initByElementTypeArray(unit_field,
+ spatial_dimension,
+ spatial_dimension,
+ _not_ghost);
+
+ ByElementTypeReal::type_iterator it = unit_field.firstType(spatial_dimension,
+ _not_ghost,
+ _ek_regular);
+ ByElementTypeReal::type_iterator end = unit_field.lastType(spatial_dimension,
+ _not_ghost,
+ _ek_regular);
+
+ for (; it != end; ++it) {
+ ElementType type = *it;
+ Array<Real> & field_array = unit_field(type);
+ UInt nb_element = mesh.getNbElement(type);
+ UInt nb_quad_per_element = model.getFEM().getNbQuadraturePoints(type);
+
+ field_array.resize(nb_element * nb_quad_per_element);
+ field_array.set(1.);
+ }
+
+ integrateFieldOnFragments(unit_field, mass);
+
+ AKANTU_DEBUG_OUT();
+}
+
+/* -------------------------------------------------------------------------- */
+void FragmentManager::computeCenterOfMass() {
+ AKANTU_DEBUG_IN();
+
+ /// integrate position multiplied by density
+ integrateFieldOnFragments(quad_coordinates, mass_center);
+
+ /// divide it by the fragments' mass
+ Real * mass_storage = mass.storage();
+ Real * mass_center_storage = mass_center.storage();
+
+ UInt total_components = mass_center.getSize() * mass_center.getNbComponent();
+
+ for (UInt i = 0; i < total_components; ++i)
+ mass_center_storage[i] /= mass_storage[i];
+
+ AKANTU_DEBUG_OUT();
+}
+
+/* -------------------------------------------------------------------------- */
+void FragmentManager::computeVelocity() {
+ AKANTU_DEBUG_IN();
+
+ UInt spatial_dimension = model.getSpatialDimension();
+
+ /// compute velocity per quadrature point
+ ByElementTypeReal velocity_field("velocity_field", id);
+
+ mesh.initByElementTypeArray(velocity_field,
+ spatial_dimension,
+ spatial_dimension,
+ _not_ghost);
+
+ model.getFEM().interpolateOnQuadraturePoints(model.getVelocity(),
+ velocity_field);
+
+ /// integrate on fragments
+ integrateFieldOnFragments(velocity_field, velocity);
+
+ /// divide it by the fragments' mass
+ Real * mass_storage = mass.storage();
+ Real * velocity_storage = velocity.storage();
+
+ UInt total_components = velocity.getSize() * velocity.getNbComponent();
+
+ for (UInt i = 0; i < total_components; ++i)
+ velocity_storage[i] /= mass_storage[i];
+
+ AKANTU_DEBUG_OUT();
+}
+
+/* -------------------------------------------------------------------------- */
+/**
+ * Given the distance @f$ \mathbf{r} @f$ between a quadrature point
+ * and its center of mass, the moment of inertia is computed as \f[
+ * I_\mathrm{CM} = \mathrm{tr}(\mathbf{r}\mathbf{r}^\mathrm{T})
+ * \mathbf{I} - \mathbf{r}\mathbf{r}^\mathrm{T} \f] for more
+ * information check Wikipedia
+ * (http://en.wikipedia.org/wiki/Moment_of_inertia#Identities_for_a_skew-symmetric_matrix)
+ *
+ */
+
+void FragmentManager::computeInertiaMoments() {
+ AKANTU_DEBUG_IN();
+
+ UInt spatial_dimension = model.getSpatialDimension();
+
+ computeCenterOfMass();
+
+ /// compute local coordinates products with respect to the center of match
+ ByElementTypeReal moments_coords("moments_coords", id);
+
+ mesh.initByElementTypeArray(moments_coords,
+ spatial_dimension * spatial_dimension,
+ spatial_dimension,
+ _not_ghost);
+
+ /// resize the by element type
+ ByElementTypeReal::type_iterator it = moments_coords.firstType(spatial_dimension,
+ _not_ghost,
+ _ek_regular);
+ ByElementTypeReal::type_iterator end = moments_coords.lastType(spatial_dimension,
+ _not_ghost,
+ _ek_regular);
+
+ for (; it != end; ++it) {
+ ElementType type = *it;
+ Array<Real> & field_array = moments_coords(type);
+ UInt nb_element = mesh.getNbElement(type);
+ UInt nb_quad_per_element = model.getFEM().getNbQuadraturePoints(type);
+
+ field_array.resize(nb_element * nb_quad_per_element);
+ }
+
+
+ /// compute coordinates
+ Array<Real>::const_iterator< Vector<Real> > mass_center_it
+ = mass_center.begin(spatial_dimension);
+
+ /// loop over fragments
+ for(const_element_group_iterator it(element_group_begin());
+ it != element_group_end(); ++it, ++mass_center_it) {
+
+ const ByElementTypeUInt & el_list = it->second->getElements();
+
+ ByElementTypeUInt::type_iterator type_it = el_list.firstType(spatial_dimension,
+ _not_ghost,
+ _ek_regular);
+ ByElementTypeUInt::type_iterator type_end = el_list.lastType(spatial_dimension,
+ _not_ghost,
+ _ek_regular);
+
+ /// loop over elements of the fragment
+ for (; type_it != type_end; ++type_it) {
+ ElementType type = *type_it;
+ UInt nb_quad_per_element = model.getFEM().getNbQuadraturePoints(type);
+
+ Array<Real> & moments_coords_array = moments_coords(type);
+ const Array<Real> & quad_coordinates_array = quad_coordinates(type);
+ const Array<UInt> & el_list_array = el_list(type);
+
+ Array<Real>::iterator< Matrix<Real> > moments_begin
+ = moments_coords_array.begin(spatial_dimension, spatial_dimension);
+ Array<Real>::const_iterator< Vector<Real> > quad_coordinates_begin
+ = quad_coordinates_array.begin(spatial_dimension);
+
+ Vector<Real> relative_coords(spatial_dimension);
+
+ for (UInt el = 0; el < el_list_array.getSize(); ++el) {
+ UInt global_el = el_list_array(el);
+
+ /// loop over quadrature points
+ for (UInt q = 0; q < nb_quad_per_element; ++q) {
+ UInt global_q = global_el * nb_quad_per_element + q;
+ Matrix<Real> & moments_matrix = moments_begin[global_q];
+ const Vector<Real> & quad_coord_vector = quad_coordinates_begin[global_q];
+
+ /// to understand this read the documentation written just
+ /// before this function
+ relative_coords = quad_coord_vector;
+ relative_coords -= *mass_center_it;
+
+ moments_matrix.outerProduct(relative_coords, relative_coords);
+ Real trace = moments_matrix.trace();
+ moments_matrix *= -1.;
+ moments_matrix += Matrix<Real>::eye(spatial_dimension, trace);
+ }
+ }
+ }
+ }
+
+ /// integrate moments
+ Array<Real> integrated_moments(global_nb_fragment,
+ spatial_dimension * spatial_dimension);
+
+ integrateFieldOnFragments(moments_coords, integrated_moments);
+
+ /// compute and store principal moments
+ inertia_moments.resize(global_nb_fragment);
+
+ Array<Real>::iterator< Matrix<Real> > integrated_moments_it
+ = integrated_moments.begin(spatial_dimension, spatial_dimension);
+ Array<Real>::iterator< Vector<Real> > inertia_moments_it
+ = inertia_moments.begin(spatial_dimension);
+
+ for (UInt frag = 0; frag < global_nb_fragment; ++frag,
+ ++integrated_moments_it, ++inertia_moments_it)
+ integrated_moments_it->eig(*inertia_moments_it);
+
+ if (dump_data)
+ createDumpDataArray(inertia_moments, "moments of inertia");
+
+ AKANTU_DEBUG_OUT();
+}
+
+/* -------------------------------------------------------------------------- */
+void FragmentManager::computeAllData() {
+ AKANTU_DEBUG_IN();
+
+ buildFragments();
+ computeVelocity();
+ computeInertiaMoments();
+
+ AKANTU_DEBUG_OUT();
+}
+
+/* -------------------------------------------------------------------------- */
+void FragmentManager::storeMassDensityPerQuadraturePoint() {
+ AKANTU_DEBUG_IN();
+
+ UInt spatial_dimension = model.getSpatialDimension();
+
+ Mesh::type_iterator it = mesh.firstType(spatial_dimension);
+ Mesh::type_iterator end = mesh.lastType(spatial_dimension);
+
+ for (; it != end; ++it) {
+ ElementType type = *it;
+
+ Array<Real> & mass_density_array = mass_density(type);
+
+ UInt nb_element = mesh.getNbElement(type);
+ UInt nb_quad_per_element = model.getFEM().getNbQuadraturePoints(type);
+ mass_density_array.resize(nb_element * nb_quad_per_element);
+
+ Array<UInt> & el_index_by_mat = model.getElementIndexByMaterial(type);
+
+ Real * mass_density_it = mass_density_array.storage();
+
+ /// store mass_density for each element and quadrature point
+ for (UInt el = 0; el < nb_element; ++el) {
+ Material & mat = model.getMaterial(el_index_by_mat(el, 0));
+
+ for (UInt q = 0; q < nb_quad_per_element; ++q, ++mass_density_it)
+ *mass_density_it = mat.getRho();
+ }
+ }
+
+ AKANTU_DEBUG_OUT();
+}
+
+/* -------------------------------------------------------------------------- */
+void FragmentManager::integrateFieldOnFragments(ByElementTypeReal & field,
+ Array<Real> & output) {
+ AKANTU_DEBUG_IN();
+
+ UInt spatial_dimension = model.getSpatialDimension();
+ UInt nb_component = output.getNbComponent();
+
+ /// integration part
+ output.resize(global_nb_fragment);
+ output.clear();
+
+ UInt * fragment_index_it = fragment_index.storage();
+ Array<Real>::iterator< Vector<Real> > output_begin = output.begin(nb_component);
+
+ /// loop over fragments
+ for(const_element_group_iterator it(element_group_begin());
+ it != element_group_end(); ++it, ++fragment_index_it) {
+
+ const ByElementTypeUInt & el_list = it->second->getElements();
+
+ ByElementTypeUInt::type_iterator type_it = el_list.firstType(spatial_dimension,
+ _not_ghost,
+ _ek_regular);
+ ByElementTypeUInt::type_iterator type_end = el_list.lastType(spatial_dimension,
+ _not_ghost,
+ _ek_regular);
+
+ /// loop over elements of the fragment
+ for (; type_it != type_end; ++type_it) {
+ ElementType type = *type_it;
+ UInt nb_quad_per_element = model.getFEM().getNbQuadraturePoints(type);
+
+ const Array<Real> & density_array = mass_density(type);
+ Array<Real> & field_array = field(type);
+ const Array<UInt> & elements = el_list(type);
+ UInt nb_element = elements.getSize();
+
+ /// generate array to be integrated by filtering fragment's elements
+ Array<Real> integration_array(elements.getSize() * nb_quad_per_element,
+ nb_component);
+
+ Array<Real>::iterator< Matrix<Real> > int_array_it
+ = integration_array.begin_reinterpret(nb_quad_per_element,
+ nb_component, nb_element);
+ Array<Real>::iterator< Matrix<Real> > int_array_end
+ = integration_array.end_reinterpret(nb_quad_per_element,
+ nb_component, nb_element);
+ Array<Real>::iterator< Matrix<Real> > field_array_begin
+ = field_array.begin_reinterpret(nb_quad_per_element,
+ nb_component,
+ field_array.getSize() / nb_quad_per_element);
+ Array<Real>::const_iterator< Vector<Real> > density_array_begin
+ = density_array.begin_reinterpret(nb_quad_per_element,
+ density_array.getSize() / nb_quad_per_element);
+
+ for (UInt el = 0; int_array_it != int_array_end; ++int_array_it, ++el) {
+ UInt global_el = elements(el);
+ *int_array_it = field_array_begin[global_el];
+
+ /// multiply field by density
+ const Vector<Real> & density_vector = density_array_begin[global_el];
+
+ for (UInt i = 0; i < nb_quad_per_element; ++i) {
+ for (UInt j = 0; j < nb_component; ++j) {
+ (*int_array_it)(i, j) *= density_vector(i);
+ }
+ }
+ }
+
+
+ /// integrate the field over the fragment
+ Array<Real> integrated_array(elements.getSize(), nb_component);
+ model.getFEM().integrate(integration_array,
+ integrated_array,
+ nb_component,
+ type,
+ _not_ghost,
+ elements);
+
+ /// sum over all elements and store the result
+ output_begin[*fragment_index_it]
+ += std::accumulate(integrated_array.begin(nb_component),
+ integrated_array.end(nb_component),
+ Vector<Real>(nb_component));
+ }
+ }
+
+ /// sum output over all processors
+ StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator();
+ comm.allReduce(output.storage(), global_nb_fragment * spatial_dimension, _so_sum);
+
+ AKANTU_DEBUG_OUT();
+}
+
+/* -------------------------------------------------------------------------- */
+template <typename T>
+void FragmentManager::createDumpDataArray(Array<T> & data, std::string name) {
+ AKANTU_DEBUG_IN();
+
+ if (data.getSize() == 0) return;
+
+ typedef typename Array<T>::template iterator< Vector<T> > data_iterator;
+ Mesh & mesh_not_const = const_cast<Mesh &>(mesh);
+
+ UInt spatial_dimension = mesh.getSpatialDimension();
+ UInt nb_component = data.getNbComponent();
+ UInt * fragment_index_it = fragment_index.storage();
+
+ data_iterator data_begin = data.begin(nb_component);
+
+ /// loop over fragments
+ for(const_element_group_iterator it(element_group_begin());
+ it != element_group_end(); ++it, ++fragment_index_it) {
+
+ const ElementGroup & fragment = *(it->second);
+
+ /// loop over cluster types
+ ElementGroup::type_iterator type_it = fragment.firstType(spatial_dimension);
+ ElementGroup::type_iterator type_end = fragment.lastType(spatial_dimension);
+
+ for(; type_it != type_end; ++type_it) {
+ ElementType type = *type_it;
+
+ /// init mesh data
+ Array<T> * mesh_data
+ = mesh_not_const.getDataPointer<T>(name, type, _not_ghost, nb_component);
+
+ data_iterator mesh_data_begin = mesh_data->begin(nb_component);
+
+ ElementGroup::const_element_iterator el_it = fragment.element_begin(type);
+ ElementGroup::const_element_iterator el_end = fragment.element_end(type);
+
+ /// fill mesh data with cluster index
+ for (; el_it != el_end; ++el_it)
+ mesh_data_begin[*el_it] = data_begin[*fragment_index_it];
+ }
+ }
+
+ AKANTU_DEBUG_OUT();
+}
+
+
+__END_AKANTU__
diff --git a/src/model/solid_mechanics/fragment_manager.hh b/src/model/solid_mechanics/fragment_manager.hh
new file mode 100644
index 000000000..757d69ad5
--- /dev/null
+++ b/src/model/solid_mechanics/fragment_manager.hh
@@ -0,0 +1,152 @@
+/**
+ * @file fragment_manager.hh
+ * @author Marco Vocialta <marco.vocialta@epfl.ch>
+ * @date Mon Jan 20 14:26:42 2014
+ *
+ * @brief Group manager to handle fragments
+ *
+ * @section LICENSE
+ *
+ * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
+ * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
+ *
+ * Akantu is free software: you can redistribute it and/or modify it under the
+ * terms of the GNU Lesser General Public License as published by the Free
+ * Software Foundation, either version 3 of the License, or (at your option) any
+ * later version.
+ *
+ * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
+ * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+ * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
+ * details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with Akantu. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+/* -------------------------------------------------------------------------- */
+
+#ifndef __AKANTU_FRAGMENT_MANAGER_HH__
+#define __AKANTU_FRAGMENT_MANAGER_HH__
+
+#include "group_manager.hh"
+#include "solid_mechanics_model_cohesive.hh"
+
+__BEGIN_AKANTU__
+
+/* -------------------------------------------------------------------------- */
+
+class FragmentManager : public GroupManager {
+ /* ------------------------------------------------------------------------ */
+ /* Constructors/Destructors */
+ /* ------------------------------------------------------------------------ */
+public:
+
+ FragmentManager(SolidMechanicsModelCohesive & model,
+ bool dump_data = true,
+ const ID & id = "fragment_manager",
+ const MemoryID & memory_id = 0);
+
+ /* ------------------------------------------------------------------------ */
+ /* Methods */
+ /* ------------------------------------------------------------------------ */
+private:
+
+ /// store mass density per quadrature point
+ void storeMassDensityPerQuadraturePoint();
+
+ /// integrate an elemental field multiplied by density on global
+ /// fragments
+ void integrateFieldOnFragments(ByElementTypeReal & field,
+ Array<Real> & output);
+
+ /// compute fragments' mass
+ void computeMass();
+
+ /// create dump data for a single array
+ template <typename T>
+ void createDumpDataArray(Array<T> & data, std::string name);
+
+public:
+
+ /// build fragment list
+ void buildFragments();
+
+ /// compute fragments' center of mass
+ void computeCenterOfMass();
+
+ /// compute fragments' velocity
+ void computeVelocity();
+
+ /// computes principal moments of inertia with respect to the center
+ /// of mass of each fragment
+ void computeInertiaMoments();
+
+ /// compute all fragments' data
+ void computeAllData();
+
+ /* ------------------------------------------------------------------------ */
+ /* Accessors */
+ /* ------------------------------------------------------------------------ */
+public:
+
+ /// get number of fragments
+ AKANTU_GET_MACRO(NbFragment, global_nb_fragment, UInt);
+
+ /// get fragments' mass
+ AKANTU_GET_MACRO(Mass, mass, const Array<Real> &);
+
+ /// get fragments' center of mass
+ AKANTU_GET_MACRO(CenterOfMass, mass_center, const Array<Real> &);
+
+ /// get fragments' velocity
+ AKANTU_GET_MACRO(Velocity, velocity, const Array<Real> &);
+
+ /// get fragments' principal moments of inertia
+ AKANTU_GET_MACRO(MomentsOfInertia, inertia_moments, const Array<Real> &);
+
+ /* ------------------------------------------------------------------------ */
+ /* Class Members */
+ /* ------------------------------------------------------------------------ */
+private:
+
+ /// local_fragment index
+ Array<UInt> fragment_index;
+
+ /// global number of fragments (parallel simulations)
+ UInt global_nb_fragment;
+
+ /// number of fragments
+ UInt nb_fragment;
+
+ /// cohesive solid mechanics model associated
+ SolidMechanicsModelCohesive & model;
+
+ /// fragments' center of mass
+ Array<Real> mass_center;
+
+ /// fragments' mass
+ Array<Real> mass;
+
+ /// fragments' velocity
+ Array<Real> velocity;
+
+ /// fragments' principal moments of inertia with respect to the
+ /// center of mass
+ Array<Real> inertia_moments;
+
+ /// quadrature points' coordinates
+ ByElementTypeReal quad_coordinates;
+
+ /// mass density per quadrature point
+ ByElementTypeReal mass_density;
+
+ /// dump data
+ bool dump_data;
+
+};
+
+__END_AKANTU__
+
+#endif /* __AKANTU_FRAGMENT_MANAGER_HH__ */
diff --git a/src/model/solid_mechanics/material.cc b/src/model/solid_mechanics/material.cc
index 6ccf33b3c..fe83e1b47 100644
--- a/src/model/solid_mechanics/material.cc
+++ b/src/model/solid_mechanics/material.cc
@@ -1,1350 +1,1351 @@
/**
* @file material.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date Tue Jul 27 18:15:37 2010
*
* @brief Implementation of the common part of the material class
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "material.hh"
#include "solid_mechanics_model.hh"
#include "sparse_matrix.hh"
#include "dof_synchronizer.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
Material::Material(SolidMechanicsModel & model, const ID & id) :
Memory(id, model.getMemoryID()),
Parsable(_st_material, id),
is_init(false),
finite_deformation(false),
inelastic_deformation(false),
name(""),
model(&model),
spatial_dimension(this->model->getSpatialDimension()),
element_filter("element_filter", id, this->memory_id),
stress("stress", *this),
strain("strain", *this),
delta_stress("delta_stress", *this),
delta_strain("delta_strain", *this),
inelas_strain("inelas_strain", *this),
piola_kirchhoff_stress("piola_kirchhoff_stress", *this),
// potential_energy_vector(false),
potential_energy("potential_energy", *this),
is_non_local(false),
use_previous_stress(false),
previous_stress("previous_stress", *this),
use_previous_strain(false),
previous_strain("previous_strain", *this),
interpolation_inverse_coordinates("interpolation inverse coordinates", *this),
interpolation_points_matrices("interpolation points matrices", *this) {
AKANTU_DEBUG_IN();
/// for each connectivity types allocate the element filer array of the material
model.getMesh().initByElementTypeArray(element_filter,
1,
spatial_dimension,
false,
_ek_regular);
registerParam("rho" , rho , 0. , _pat_parsable | _pat_modifiable, "Density");
registerParam("name" , name , std::string(), _pat_parsable | _pat_readable);
registerParam("finite_deformation" , finite_deformation , false , _pat_parsable | _pat_modifiable, "Is finite deformation");
registerParam("inelastic_deformation", inelastic_deformation, false , _pat_parsable | _pat_modifiable, "Is inelastic deformation");
/// allocate strain stress for local elements
strain.initialize(spatial_dimension * spatial_dimension);
stress.initialize(spatial_dimension * spatial_dimension);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
Material::~Material() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void Material::initMaterial() {
AKANTU_DEBUG_IN();
if(finite_deformation) {
this->delta_strain.initialize(spatial_dimension * spatial_dimension);
this->delta_stress.initialize(spatial_dimension * spatial_dimension);
this->piola_kirchhoff_stress.initialize(spatial_dimension * spatial_dimension);
}
if(inelastic_deformation) {
this->delta_strain.initialize(spatial_dimension * spatial_dimension);
this->delta_stress.initialize(spatial_dimension * spatial_dimension);
this->inelas_strain.initialize(spatial_dimension * spatial_dimension);
}
if(use_previous_stress) {
this->previous_stress.initialize(spatial_dimension * spatial_dimension);
}
if(use_previous_strain) {
this->previous_strain.initialize(spatial_dimension * spatial_dimension);
}
for (std::map<ID, InternalField<Real> *>::iterator it = internal_vectors_real.begin();
it != internal_vectors_real.end();
++it) it->second->resize();
is_init = true;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void Material::savePreviousState(const GhostType ghost_type) {
AKANTU_DEBUG_IN();
UInt spatial_dimension = model->getSpatialDimension();
Mesh::type_iterator it = model->getMesh().firstType(spatial_dimension, ghost_type);
Mesh::type_iterator last_type = model->getMesh().lastType(spatial_dimension, ghost_type);
if(finite_deformation)
for (; it != last_type; ++it) {
if(use_previous_stress) previous_stress(*it, ghost_type).copy(piola_kirchhoff_stress(*it, ghost_type));
if(use_previous_strain) previous_strain(*it, ghost_type).copy(strain(*it, ghost_type));
}
else
for (; it != last_type; ++it) {
if(use_previous_stress) previous_stress(*it, ghost_type).copy(stress(*it, ghost_type));
if(use_previous_strain) previous_strain(*it, ghost_type).copy(strain(*it, ghost_type));
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Compute the residual by assembling @f$\int_{e} \sigma_e \frac{\partial
* \varphi}{\partial X} dX @f$
*
* @param[in] displacements nodes displacements
* @param[in] ghost_type compute the residual for _ghost or _not_ghost element
*/
void Material::updateResidual(GhostType ghost_type) {
AKANTU_DEBUG_IN();
computeAllStresses(ghost_type);
assembleResidual(ghost_type);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void Material::assembleResidual(GhostType ghost_type) {
AKANTU_DEBUG_IN();
UInt spatial_dimension = model->getSpatialDimension();
if(!finite_deformation){
Array<Real> & residual = const_cast<Array<Real> &>(model->getResidual());
Mesh & mesh = model->getFEM().getMesh();
Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type);
Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type);
for(; it != last_type; ++it) {
const Array<Real> & shapes_derivatives = model->getFEM().getShapesDerivatives(*it, ghost_type);
Array<UInt> & elem_filter = element_filter(*it, ghost_type);
UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent();
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it);
UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(*it, ghost_type);
UInt nb_element = elem_filter.getSize();
/// compute @f$\sigma \frac{\partial \varphi}{\partial X}@f$ by @f$\mathbf{B}^t \mathbf{\sigma}_q@f$
Array<Real> * sigma_dphi_dx =
new Array<Real>(nb_element*nb_quadrature_points,
size_of_shapes_derivatives, "sigma_x_dphi_/_dX");
Array<Real> * shapesd_filtered =
new Array<Real>(0, size_of_shapes_derivatives, "filtered shapesd");
FEM::filterElementalData(mesh, shapes_derivatives, *shapesd_filtered,
*it, ghost_type, elem_filter);
Array<Real> & stress_vect = stress(*it, ghost_type);
Array<Real>::iterator< Matrix<Real> > sigma =
stress_vect.begin(spatial_dimension, spatial_dimension);
Array<Real>::iterator< Matrix<Real> > B =
shapesd_filtered->begin(spatial_dimension, nb_nodes_per_element);
Array<Real>::iterator< Matrix<Real> > Bt_sigma_it =
sigma_dphi_dx->begin(spatial_dimension, nb_nodes_per_element);
for (UInt q = 0; q < nb_element*nb_quadrature_points; ++q, ++sigma, ++B, ++Bt_sigma_it)
Bt_sigma_it->mul<false,false>(*sigma, *B);
delete shapesd_filtered;
/**
* compute @f$\int \sigma * \frac{\partial \varphi}{\partial X}dX@f$ by @f$ \sum_q \mathbf{B}^t
* \mathbf{\sigma}_q \overline w_q J_q@f$
*/
Array<Real> * int_sigma_dphi_dx = new Array<Real>(nb_element,
nb_nodes_per_element * spatial_dimension,
"int_sigma_x_dphi_/_dX");
model->getFEM().integrate(*sigma_dphi_dx, *int_sigma_dphi_dx,
size_of_shapes_derivatives,
*it, ghost_type,
elem_filter);
delete sigma_dphi_dx;
/// assemble
model->getFEM().assembleArray(*int_sigma_dphi_dx, residual,
model->getDOFSynchronizer().getLocalDOFEquationNumbers(),
residual.getNbComponent(),
*it, ghost_type, elem_filter, -1);
delete int_sigma_dphi_dx;
}
}
else{
switch (spatial_dimension){
case 1:
this->assembleResidual<1>(ghost_type);
break;
case 2:
this->assembleResidual<2>(ghost_type);
break;
case 3:
this->assembleResidual<3>(ghost_type);
break;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Compute the stress from the strain
*
* @param[in] current_position nodes postition + displacements
* @param[in] ghost_type compute the residual for _ghost or _not_ghost element
*/
void Material::computeAllStresses(GhostType ghost_type) {
AKANTU_DEBUG_IN();
UInt spatial_dimension = model->getSpatialDimension();
Mesh::type_iterator it = model->getFEM().getMesh().firstType(spatial_dimension, ghost_type);
Mesh::type_iterator last_type = model->getFEM().getMesh().lastType(spatial_dimension, ghost_type);
for(; it != last_type; ++it) {
Array<UInt> & elem_filter = element_filter(*it, ghost_type);
Array<Real> & strain_vect = strain(*it, ghost_type);
/// compute @f$\nabla u@f$
model->getFEM().gradientOnQuadraturePoints(model->getDisplacement(), strain_vect,
spatial_dimension,
*it, ghost_type, elem_filter);
if(finite_deformation || inelastic_deformation){ /// compute @f$\nabla \delta u@f$
Array<Real> & delta_strain_vect = delta_strain(*it, ghost_type);
- model->getFEM().gradientOnQuadraturePoints(model->getIncrement(), delta_strain_vect,
- spatial_dimension,
- *it, ghost_type, elem_filter);
+ model->getFEM().gradientOnQuadraturePoints(model->getIncrement(),
+ delta_strain_vect,
+ spatial_dimension,
+ *it, ghost_type, elem_filter);
}
/// compute @f$\mathbf{\sigma}_q@f$ from @f$\nabla u@f$
computeStress(*it, ghost_type);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void Material::computeAllCauchyStresses(GhostType ghost_type) {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(finite_deformation,"The Cauchy stress can only be computed if you are working in finite deformation.");
//resizeInternalArray(stress);
Mesh::type_iterator it = model->getFEM().getMesh().firstType(spatial_dimension, ghost_type);
Mesh::type_iterator last_type = model->getFEM().getMesh().lastType(spatial_dimension, ghost_type);
for(; it != last_type; ++it)
switch (spatial_dimension){
case 1:
this->computeCauchyStress<1>(*it, ghost_type);
break;
case 2:
this->computeCauchyStress<2>(*it, ghost_type);
break;
case 3:
this->computeCauchyStress<3>(*it, ghost_type);
break;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
void Material::computeCauchyStress(ElementType el_type, GhostType ghost_type) {
AKANTU_DEBUG_IN();
Array<Real>::iterator< Matrix<Real> > strain_it =
this->strain(el_type, ghost_type).begin(dim, dim);
Array<Real>::iterator< Matrix<Real> > strain_end =
this->strain(el_type, ghost_type).end(dim, dim);
Array<Real>::iterator< Matrix<Real> > piola_it =
this->piola_kirchhoff_stress(el_type, ghost_type).begin(dim, dim);
Array<Real>::iterator< Matrix<Real> > stress_it =
this->stress(el_type, ghost_type).begin(dim, dim);
Matrix<Real> F_tensor(dim, dim);
for (; strain_it != strain_end; ++strain_it, ++piola_it, ++stress_it) {
Matrix<Real> & grad_u = *strain_it;
Matrix<Real> & piola = *piola_it;
Matrix<Real> & sigma = *stress_it;
gradUToF<dim > (grad_u, F_tensor);
computeCauchyStressOnQuad<dim >(F_tensor, piola, sigma);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void Material::setToSteadyState(GhostType ghost_type) {
AKANTU_DEBUG_IN();
const Array<Real> & displacement = model->getDisplacement();
//resizeInternalArray(strain);
UInt spatial_dimension = model->getSpatialDimension();
Mesh::type_iterator it = model->getFEM().getMesh().firstType(spatial_dimension, ghost_type);
Mesh::type_iterator last_type = model->getFEM().getMesh().lastType(spatial_dimension, ghost_type);
for(; it != last_type; ++it) {
Array<UInt> & elem_filter = element_filter(*it, ghost_type);
Array<Real> & strain_vect = strain(*it, ghost_type);
/// compute @f$\nabla u@f$
model->getFEM().gradientOnQuadraturePoints(displacement, strain_vect,
spatial_dimension,
*it, ghost_type, elem_filter);
setToSteadyState(*it, ghost_type);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Compute the stiffness matrix by assembling @f$\int_{\omega} B^t \times D
* \times B d\omega @f$
*
* @param[in] current_position nodes postition + displacements
* @param[in] ghost_type compute the residual for _ghost or _not_ghost element
*/
void Material::assembleStiffnessMatrix(GhostType ghost_type) {
AKANTU_DEBUG_IN();
UInt spatial_dimension = model->getSpatialDimension();
Mesh & mesh = model->getFEM().getMesh();
Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type);
Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type);
for(; it != last_type; ++it) {
if(finite_deformation){
switch (spatial_dimension) {
case 1:
{
assembleStiffnessMatrixNL < 1 > (*it, ghost_type);
assembleStiffnessMatrixL2 < 1 > (*it, ghost_type);
break;
}
case 2:
{
assembleStiffnessMatrixNL < 2 > (*it, ghost_type);
assembleStiffnessMatrixL2 < 2 > (*it, ghost_type);
break;
}
case 3:
{
assembleStiffnessMatrixNL < 3 > (*it, ghost_type);
assembleStiffnessMatrixL2 < 3 > (*it, ghost_type);
break;
}
}
} else {
switch(spatial_dimension) {
case 1: { assembleStiffnessMatrix<1>(*it, ghost_type); break; }
case 2: { assembleStiffnessMatrix<2>(*it, ghost_type); break; }
case 3: { assembleStiffnessMatrix<3>(*it, ghost_type); break; }
}
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
void Material::assembleStiffnessMatrix(const ElementType & type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
SparseMatrix & K = const_cast<SparseMatrix &>(model->getStiffnessMatrix());
const Array<Real> & shapes_derivatives = model->getFEM().getShapesDerivatives(type,ghost_type);
Array<UInt> & elem_filter = element_filter(type, ghost_type);
Array<Real> & strain_vect = strain(type, ghost_type);
UInt nb_element = elem_filter.getSize();
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(type, ghost_type);
strain_vect.resize(nb_quadrature_points * nb_element);
model->getFEM().gradientOnQuadraturePoints(model->getDisplacement(), strain_vect,
dim, type, ghost_type, elem_filter);
if(inelastic_deformation){ /// compute @f$\nabla \delta u@f$
Array<Real> & delta_strain_vect = delta_strain(type, ghost_type);
delta_strain_vect.resize(nb_quadrature_points * nb_element);
model->getFEM().gradientOnQuadraturePoints(model->getIncrement(), delta_strain_vect,
dim,
type, ghost_type, elem_filter);
}
UInt tangent_size = getTangentStiffnessVoigtSize(dim);
Array<Real> * tangent_stiffness_matrix =
new Array<Real>(nb_element*nb_quadrature_points, tangent_size * tangent_size,
"tangent_stiffness_matrix");
tangent_stiffness_matrix->clear();
computeTangentModuli(type, *tangent_stiffness_matrix, ghost_type);
Array<Real> * shapesd_filtered =
new Array<Real>(0, dim * nb_nodes_per_element, "filtered shapesd");
FEM::filterElementalData(model->getFEM().getMesh(), shapes_derivatives, *shapesd_filtered,
type, ghost_type, elem_filter);
/// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$
UInt bt_d_b_size = dim * nb_nodes_per_element;
Array<Real> * bt_d_b = new Array<Real>(nb_element * nb_quadrature_points,
bt_d_b_size * bt_d_b_size,
"B^t*D*B");
Matrix<Real> B(tangent_size, dim * nb_nodes_per_element);
Matrix<Real> Bt_D(dim * nb_nodes_per_element, tangent_size);
Array<Real>::iterator< Matrix<Real> > shapes_derivatives_filtered_it =
shapesd_filtered->begin(dim, nb_nodes_per_element);
Array<Real>::iterator< Matrix<Real> > Bt_D_B_it = bt_d_b->begin(dim*nb_nodes_per_element,
dim*nb_nodes_per_element);
Array<Real>::iterator< Matrix<Real> > D_it = tangent_stiffness_matrix->begin(tangent_size,
tangent_size);
Array<Real>::iterator< Matrix<Real> > D_end = tangent_stiffness_matrix->end (tangent_size,
tangent_size);
for(; D_it != D_end; ++D_it, ++Bt_D_B_it, ++shapes_derivatives_filtered_it) {
Matrix<Real> & D = *D_it;
Matrix<Real> & Bt_D_B = *Bt_D_B_it;
VoigtHelper<dim>::transferBMatrixToSymVoigtBMatrix(*shapes_derivatives_filtered_it, B, nb_nodes_per_element);
Bt_D.mul<true, false>(B, D);
Bt_D_B.mul<false, false>(Bt_D, B);
}
delete tangent_stiffness_matrix;
delete shapesd_filtered;
/// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$
Array<Real> * K_e = new Array<Real>(nb_element,
bt_d_b_size * bt_d_b_size,
"K_e");
model->getFEM().integrate(*bt_d_b, *K_e,
bt_d_b_size * bt_d_b_size,
type, ghost_type,
elem_filter);
delete bt_d_b;
model->getFEM().assembleMatrix(*K_e, K, spatial_dimension, type, ghost_type, elem_filter);
delete K_e;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
void Material::assembleStiffnessMatrixNL(const ElementType & type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
SparseMatrix & K = const_cast<SparseMatrix &> (model->getStiffnessMatrix());
const Array<Real> & shapes_derivatives = model->getFEM().getShapesDerivatives(type, ghost_type);
Array<UInt> & elem_filter = element_filter(type, ghost_type);
//Array<Real> & strain_vect = delta_strain(type, ghost_type);
UInt nb_element = elem_filter.getSize();
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(type, ghost_type);
//strain_vect.resize(nb_quadrature_points * nb_element);
// model->getFEM().gradientOnQuadraturePoints(model->getIncrement(), strain_vect,
// dim, type, ghost_type, &elem_filter);
Array<Real> * shapes_derivatives_filtered = new Array<Real > (nb_element * nb_quadrature_points,
dim * nb_nodes_per_element,
"shapes derivatives filtered");
Array<Real>::const_iterator< Matrix<Real> > shapes_derivatives_it = shapes_derivatives.begin(spatial_dimension,
nb_nodes_per_element);
Array<Real>::iterator< Matrix<Real> > shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(spatial_dimension,
nb_nodes_per_element);
UInt * elem_filter_val = elem_filter.storage();
for (UInt e = 0; e < nb_element; ++e, ++elem_filter_val)
for (UInt q = 0; q < nb_quadrature_points; ++q, ++shapes_derivatives_filtered_it)
*shapes_derivatives_filtered_it = shapes_derivatives_it[*elem_filter_val * nb_quadrature_points + q];
/// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$
UInt bt_s_b_size = dim * nb_nodes_per_element;
Array<Real> * bt_s_b = new Array<Real > (nb_element * nb_quadrature_points,
bt_s_b_size * bt_s_b_size,
"B^t*D*B");
UInt piola_matrix_size = getCauchyStressMatrixSize(dim);
Matrix<Real> B(piola_matrix_size, bt_s_b_size);
Matrix<Real> Bt_S(bt_s_b_size, piola_matrix_size);
Matrix<Real> S(piola_matrix_size, piola_matrix_size);
shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(spatial_dimension, nb_nodes_per_element);
Array<Real>::iterator< Matrix<Real> > Bt_S_B_it = bt_s_b->begin(bt_s_b_size,
bt_s_b_size);
Array<Real>::iterator< Matrix<Real> > Bt_S_B_end = bt_s_b->end(bt_s_b_size,
bt_s_b_size);
Array<Real>::iterator< Matrix<Real> > piola_it = piola_kirchhoff_stress(type, ghost_type).begin(dim, dim);
for (; Bt_S_B_it != Bt_S_B_end; ++Bt_S_B_it, ++shapes_derivatives_filtered_it, ++piola_it) {
Matrix<Real> & Bt_S_B = *Bt_S_B_it;
Matrix<Real> & Piola_kirchhoff_matrix = *piola_it;
SetCauchyStressMatrix< dim >(Piola_kirchhoff_matrix, S);
VoigtHelper<dim>::transferBMatrixToBNL(*shapes_derivatives_filtered_it, B, nb_nodes_per_element);
Bt_S.mul < true, false > (B, S);
Bt_S_B.mul < false, false > (Bt_S, B);
}
delete shapes_derivatives_filtered;
/// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$
Array<Real> * K_e = new Array<Real > (nb_element,
bt_s_b_size * bt_s_b_size,
"K_e");
model->getFEM().integrate(*bt_s_b, *K_e,
bt_s_b_size * bt_s_b_size,
type, ghost_type,
elem_filter);
delete bt_s_b;
model->getFEM().assembleMatrix(*K_e, K, spatial_dimension, type, ghost_type, elem_filter);
delete K_e;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
void Material::assembleStiffnessMatrixL2(const ElementType & type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
SparseMatrix & K = const_cast<SparseMatrix &> (model->getStiffnessMatrix());
const Array<Real> & shapes_derivatives = model->getFEM().getShapesDerivatives(type, ghost_type);
Array<UInt> & elem_filter = element_filter(type, ghost_type);
Array<Real> & strain_vect = strain(type, ghost_type);
UInt nb_element = elem_filter.getSize();
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(type, ghost_type);
strain_vect.resize(nb_quadrature_points * nb_element);
model->getFEM().gradientOnQuadraturePoints(model->getDisplacement(), strain_vect,
dim, type, ghost_type, elem_filter);
UInt tangent_size = getTangentStiffnessVoigtSize(dim);
Array<Real> * tangent_stiffness_matrix =
new Array<Real > (nb_element*nb_quadrature_points, tangent_size * tangent_size,
"tangent_stiffness_matrix");
tangent_stiffness_matrix->clear();
computeTangentModuli(type, *tangent_stiffness_matrix, ghost_type);
Array<Real> * shapes_derivatives_filtered = new Array<Real > (nb_element * nb_quadrature_points,
dim * nb_nodes_per_element,
"shapes derivatives filtered");
Array<Real>::const_iterator< Matrix<Real> > shapes_derivatives_it = shapes_derivatives.begin(spatial_dimension,
nb_nodes_per_element);
Array<Real>::iterator< Matrix<Real> > shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(spatial_dimension,
nb_nodes_per_element);
UInt * elem_filter_val = elem_filter.storage();
for (UInt e = 0; e < nb_element; ++e, ++elem_filter_val)
for (UInt q = 0; q < nb_quadrature_points; ++q, ++shapes_derivatives_filtered_it)
*shapes_derivatives_filtered_it = shapes_derivatives_it[*elem_filter_val * nb_quadrature_points + q];
/// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$
UInt bt_d_b_size = dim * nb_nodes_per_element;
Array<Real> * bt_d_b = new Array<Real > (nb_element * nb_quadrature_points,
bt_d_b_size * bt_d_b_size,
"B^t*D*B");
Matrix<Real> B(tangent_size, dim * nb_nodes_per_element);
Matrix<Real> B2(tangent_size, dim * nb_nodes_per_element);
Matrix<Real> Bt_D(dim * nb_nodes_per_element, tangent_size);
shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(spatial_dimension, nb_nodes_per_element);
Array<Real>::iterator< Matrix<Real> > Bt_D_B_it = bt_d_b->begin(dim*nb_nodes_per_element,
dim * nb_nodes_per_element);
Array<Real>::iterator< Matrix<Real> > grad_u_it = strain_vect.begin(dim, dim);
Array<Real>::iterator< Matrix<Real> > D_it = tangent_stiffness_matrix->begin(tangent_size,
tangent_size);
Array<Real>::iterator< Matrix<Real> > D_end = tangent_stiffness_matrix->end(tangent_size,
tangent_size);
for (; D_it != D_end; ++D_it, ++Bt_D_B_it, ++shapes_derivatives_filtered_it, ++grad_u_it) {
Matrix<Real> & grad_u = *grad_u_it;
Matrix<Real> & D = *D_it;
Matrix<Real> & Bt_D_B = *Bt_D_B_it;
//transferBMatrixToBL1<dim > (*shapes_derivatives_filtered_it, B, nb_nodes_per_element);
VoigtHelper<dim>::transferBMatrixToSymVoigtBMatrix(*shapes_derivatives_filtered_it, B, nb_nodes_per_element);
VoigtHelper<dim>::transferBMatrixToBL2(*shapes_derivatives_filtered_it, grad_u, B2, nb_nodes_per_element);
B += B2;
Bt_D.mul < true, false > (B, D);
Bt_D_B.mul < false, false > (Bt_D, B);
}
delete tangent_stiffness_matrix;
delete shapes_derivatives_filtered;
/// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$
Array<Real> * K_e = new Array<Real > (nb_element,
bt_d_b_size * bt_d_b_size,
"K_e");
model->getFEM().integrate(*bt_d_b, *K_e,
bt_d_b_size * bt_d_b_size,
type, ghost_type,
elem_filter);
delete bt_d_b;
model->getFEM().assembleMatrix(*K_e, K, spatial_dimension, type, ghost_type, elem_filter);
delete K_e;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
void Material::assembleResidual(GhostType ghost_type){
AKANTU_DEBUG_IN();
Array<Real> & residual = const_cast<Array<Real> &> (model->getResidual());
Mesh & mesh = model->getFEM().getMesh();
Mesh::type_iterator it = mesh.firstType(dim, ghost_type);
Mesh::type_iterator last_type = mesh.lastType(dim, ghost_type);
for (; it != last_type; ++it) {
const Array<Real> & shapes_derivatives = model->getFEM().getShapesDerivatives(*it, ghost_type);
Array<UInt> & elem_filter = element_filter(*it, ghost_type);
UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent();
UInt nb_element = elem_filter.getSize();
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it);
UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(*it, ghost_type);
Array<Real> * shapesd_filtered =
new Array<Real>(0, size_of_shapes_derivatives, "filtered shapesd");
FEM::filterElementalData(mesh, shapes_derivatives, *shapesd_filtered,
*it, ghost_type, elem_filter);
Array<Real>::iterator< Matrix<Real> > shapes_derivatives_filtered_it = shapesd_filtered->begin(dim,
nb_nodes_per_element);
//Set stress vectors
UInt stress_size = getTangentStiffnessVoigtSize(dim);
//Set matrices B and BNL*
UInt bt_s_size = dim * nb_nodes_per_element;
Array<Real> * bt_s = new Array<Real > (nb_element * nb_quadrature_points,
bt_s_size,
"B^t*S");
Array<Real>::iterator< Matrix<Real> > grad_u_it = this->strain(*it, ghost_type).begin(dim, dim);
Array<Real>::iterator< Matrix<Real> > grad_u_end = this->strain(*it, ghost_type).end(dim, dim);
Array<Real>::iterator< Matrix<Real> > stress_it = this->piola_kirchhoff_stress(*it, ghost_type).begin(dim, dim);
shapes_derivatives_filtered_it = shapesd_filtered->begin(dim, nb_nodes_per_element);
Array<Real>::iterator< Matrix<Real> > bt_s_it = bt_s->begin(bt_s_size, 1);
Matrix<Real> S_vect(stress_size, 1);
Matrix<Real> B_tensor(stress_size, bt_s_size);
Matrix<Real> B2_tensor(stress_size, bt_s_size);
for (; grad_u_it != grad_u_end; ++grad_u_it, ++stress_it, ++shapes_derivatives_filtered_it, ++bt_s_it) {
Matrix<Real> & grad_u = *grad_u_it;
Matrix<Real> & r_it = *bt_s_it;
Matrix<Real> & S_it = *stress_it;
SetCauchyStressArray <dim> (S_it, S_vect);
VoigtHelper<dim>::transferBMatrixToSymVoigtBMatrix(*shapes_derivatives_filtered_it, B_tensor, nb_nodes_per_element);
VoigtHelper<dim>::transferBMatrixToBL2(*shapes_derivatives_filtered_it, grad_u, B2_tensor, nb_nodes_per_element);
B_tensor += B2_tensor;
r_it.mul < true, false > (B_tensor, S_vect);
}
delete shapesd_filtered;
/// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$
Array<Real> * r_e = new Array<Real > (nb_element,
bt_s_size, "r_e");
model->getFEM().integrate(*bt_s, *r_e,
bt_s_size,
*it, ghost_type,
elem_filter);
delete bt_s;
model->getFEM().assembleArray(*r_e, residual,
model->getDOFSynchronizer().getLocalDOFEquationNumbers(),
residual.getNbComponent(),
*it, ghost_type, elem_filter, -1);
delete r_e;
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void Material::computeAllStressesFromTangentModuli(GhostType ghost_type) {
AKANTU_DEBUG_IN();
UInt spatial_dimension = model->getSpatialDimension();
Mesh & mesh = model->getFEM().getMesh();
Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type);
Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type);
for(; it != last_type; ++it) {
switch(spatial_dimension) {
case 1: { computeAllStressesFromTangentModuli<1>(*it, ghost_type); break; }
case 2: { computeAllStressesFromTangentModuli<2>(*it, ghost_type); break; }
case 3: { computeAllStressesFromTangentModuli<3>(*it, ghost_type); break; }
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template<UInt dim>
void Material::computeAllStressesFromTangentModuli(const ElementType & type,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
const Array<Real> & shapes_derivatives = model->getFEM().getShapesDerivatives(type, ghost_type);
Array<UInt> & elem_filter = element_filter(type, ghost_type);
Array<Real> & strain_vect = strain(type, ghost_type);
UInt nb_element = elem_filter.getSize();
UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type);
UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(type, ghost_type);
strain_vect.resize(nb_quadrature_points * nb_element);
Array<Real> & disp = model->getDisplacement();
model->getFEM().gradientOnQuadraturePoints(disp, strain_vect,
dim, type, ghost_type, elem_filter);
UInt tangent_moduli_size = getTangentStiffnessVoigtSize(dim);
Array<Real> * tangent_moduli_tensors =
new Array<Real>(nb_element*nb_quadrature_points, tangent_moduli_size * tangent_moduli_size,
"tangent_moduli_tensors");
tangent_moduli_tensors->clear();
computeTangentModuli(type, *tangent_moduli_tensors, ghost_type);
Array<Real> * shapesd_filtered =
new Array<Real>(0, dim* nb_nodes_per_element, "filtered shapesd");
FEM::filterElementalData(model->getFEM().getMesh(), shapes_derivatives, *shapesd_filtered,
type, ghost_type, elem_filter);
Array<Real> filtered_u(nb_element, nb_nodes_per_element * spatial_dimension);
FEM::extractNodalToElementField(model->getFEM().getMesh(), disp, filtered_u,
type, ghost_type, elem_filter);
/// compute @f$\mathbf{D} \mathbf{B} \mathbf{u}@f$
Array<Real>::iterator< Matrix<Real> > shapes_derivatives_filtered_it =
shapesd_filtered->begin(dim, nb_nodes_per_element);
Array<Real>::iterator< Matrix<Real> > D_it = tangent_moduli_tensors->begin(tangent_moduli_size,
tangent_moduli_size);
Array<Real>::iterator< Matrix<Real> > sigma_it = stress(type, ghost_type).begin(spatial_dimension,
spatial_dimension);
Array<Real>::iterator< Vector<Real> > u_it = filtered_u.begin(spatial_dimension * nb_nodes_per_element);
Matrix<Real> B(tangent_moduli_size, spatial_dimension * nb_nodes_per_element);
Vector<Real> Bu(tangent_moduli_size);
Vector<Real> DBu(tangent_moduli_size);
for (UInt e = 0; e < nb_element; ++e, ++u_it) {
for (UInt q = 0; q < nb_quadrature_points; ++q, ++D_it, ++shapes_derivatives_filtered_it, ++sigma_it) {
Vector<Real> & u = *u_it;
Matrix<Real> & sigma = *sigma_it;
Matrix<Real> & D = *D_it;
VoigtHelper<dim>::transferBMatrixToSymVoigtBMatrix(*shapes_derivatives_filtered_it, B, nb_nodes_per_element);
Bu.mul<false>(B, u);
DBu.mul<false>(D, Bu);
// Voigt notation to full symmetric tensor
for (UInt i = 0; i < dim; ++i) sigma(i, i) = DBu(i);
if(dim == 2) {
sigma(0,1) = sigma(1,0) = DBu(2);
} else if(dim == 3) {
sigma(1,2) = sigma(2,1) = DBu(3);
sigma(0,2) = sigma(2,0) = DBu(4);
sigma(0,1) = sigma(1,0) = DBu(5);
}
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void Material::computePotentialEnergy(ElementType el_type, GhostType ghost_type) {
AKANTU_DEBUG_IN();
if(ghost_type != _not_ghost) return;
Real * epot = potential_energy(el_type, ghost_type).storage();
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type);
computePotentialEnergyOnQuad(grad_u, sigma, *epot);
epot++;
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void Material::computePotentialEnergyByElement() {
AKANTU_DEBUG_IN();
Mesh & mesh = model->getFEM().getMesh();
Mesh::type_iterator it = mesh.firstType(spatial_dimension);
Mesh::type_iterator last_type = mesh.lastType(spatial_dimension);
for(; it != last_type; ++it) {
if(!potential_energy.exists(*it, _not_ghost)) {
UInt nb_element = element_filter(*it, _not_ghost).getSize();
UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(*it, _not_ghost);
potential_energy.alloc(nb_element * nb_quadrature_points, 1,
*it, _not_ghost);
}
computePotentialEnergy(*it);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void Material::computePotentialEnergyByElement(ElementType type, UInt index,
Vector<Real> & epot_on_quad_points){
Array<Real>::iterator< Matrix<Real> > strain_it =
this->strain(type).begin(spatial_dimension,
spatial_dimension);
Array<Real>::iterator< Matrix<Real> > strain_end =
this->strain(type).begin(spatial_dimension,
spatial_dimension);
Array<Real>::iterator< Matrix<Real> > stress_it =
this->stress(type).begin(spatial_dimension,
spatial_dimension);
UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(type);
strain_it += index*nb_quadrature_points;
strain_end += (index+1)*nb_quadrature_points;
stress_it += index*nb_quadrature_points;
Real * epot_quad = epot_on_quad_points.storage();
for(;strain_it != strain_end; ++strain_it, ++stress_it, ++epot_quad) {
Matrix<Real> & grad_u = *strain_it;
Matrix<Real> & sigma = *stress_it;
computePotentialEnergyOnQuad(grad_u, sigma, *epot_quad);
}
}
/* -------------------------------------------------------------------------- */
Real Material::getPotentialEnergy() {
AKANTU_DEBUG_IN();
Real epot = 0.;
computePotentialEnergyByElement();
/// integrate the potential energy for each type of elements
Mesh & mesh = model->getFEM().getMesh();
Mesh::type_iterator it = mesh.firstType(spatial_dimension);
Mesh::type_iterator last_type = mesh.lastType(spatial_dimension);
for(; it != last_type; ++it) {
epot += model->getFEM().integrate(potential_energy(*it, _not_ghost), *it,
_not_ghost, element_filter(*it, _not_ghost));
}
AKANTU_DEBUG_OUT();
return epot;
}
/* -------------------------------------------------------------------------- */
Real Material::getPotentialEnergy(ElementType & type, UInt index) {
AKANTU_DEBUG_IN();
Real epot = 0.;
Vector<Real> epot_on_quad_points(model->getFEM().getNbQuadraturePoints(type));
computePotentialEnergyByElement(type,index,epot_on_quad_points);
epot = model->getFEM().integrate(epot_on_quad_points, type, element_filter(type)(index));
AKANTU_DEBUG_OUT();
return epot;
}
/* -------------------------------------------------------------------------- */
Real Material::getEnergy(std::string type) {
AKANTU_DEBUG_IN();
if(type == "potential") return getPotentialEnergy();
AKANTU_DEBUG_OUT();
return 0.;
}
/* -------------------------------------------------------------------------- */
Real Material::getEnergy(std::string energy_id, ElementType type, UInt index) {
AKANTU_DEBUG_IN();
if(energy_id == "potential") return getPotentialEnergy(type, index);
AKANTU_DEBUG_OUT();
return 0.;
}
/* -------------------------------------------------------------------------- */
void Material::computeQuadraturePointsCoordinates(ByElementTypeReal & quadrature_points_coordinates, const GhostType & ghost_type) const {
AKANTU_DEBUG_IN();
const Mesh & mesh = model->getFEM().getMesh();
Array<Real> nodes_coordinates(mesh.getNodes(), true);
nodes_coordinates += model->getDisplacement();
Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type);
Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type);
for(; it != last_type; ++it) {
const Array<UInt> & elem_filter = element_filter(*it, ghost_type);
UInt nb_element = elem_filter.getSize();
UInt nb_tot_quad = model->getFEM().getNbQuadraturePoints(*it, ghost_type) * nb_element;
Array<Real> & quads = quadrature_points_coordinates(*it, ghost_type);
quads.resize(nb_tot_quad);
model->getFEM().interpolateOnQuadraturePoints(nodes_coordinates,
quads, spatial_dimension,
*it, ghost_type, elem_filter);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void Material::initElementalFieldInterpolation(const ByElementTypeReal & interpolation_points_coordinates) {
AKANTU_DEBUG_IN();
const Mesh & mesh = model->getFEM().getMesh();
ByElementTypeArray<Real> quadrature_points_coordinates("quadrature_points_coordinates_for_interpolation", getID());
mesh.initByElementTypeArray(quadrature_points_coordinates, spatial_dimension, spatial_dimension);
for (ghost_type_t::iterator gt = ghost_type_t::begin();
gt != ghost_type_t::end(); ++gt) {
GhostType ghost_type = *gt;
computeQuadraturePointsCoordinates(quadrature_points_coordinates, ghost_type);
Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type);
Mesh::type_iterator last = mesh.lastType(spatial_dimension, ghost_type);
for (; it != last; ++it) {
ElementType type = *it;
UInt nb_element = mesh.getNbElement(type, ghost_type);
if (nb_element == 0) continue;
const Array<Real> & interp_points_coord = interpolation_points_coordinates(type, ghost_type);
UInt nb_interpolation_points_per_elem = interp_points_coord.getSize() / nb_element;
AKANTU_DEBUG_ASSERT(interp_points_coord.getSize() % nb_element == 0,
"Number of interpolation points is wrong");
#define AKANTU_INIT_INTERPOLATE_ELEMENTAL_FIELD(type) \
initElementalFieldInterpolation<type>(quadrature_points_coordinates(type, ghost_type), \
interp_points_coord, \
nb_interpolation_points_per_elem, \
ghost_type) \
AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(AKANTU_INIT_INTERPOLATE_ELEMENTAL_FIELD);
#undef AKANTU_INIT_INTERPOLATE_ELEMENTAL_FIELD
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <ElementType type>
void Material::initElementalFieldInterpolation(const Array<Real> & quad_coordinates,
const Array<Real> & interpolation_points_coordinates,
const UInt nb_interpolation_points_per_elem,
const GhostType ghost_type) {
AKANTU_DEBUG_IN();
UInt size_inverse_coords =
getSizeElementalFieldInterpolationCoodinates<type>(ghost_type);
Array<UInt> & elem_fil = element_filter(type, ghost_type);
UInt nb_element = elem_fil.getSize();
UInt nb_quad_per_element = model->getFEM().getNbQuadraturePoints(type, ghost_type);
if(!interpolation_inverse_coordinates.exists(type, ghost_type))
interpolation_inverse_coordinates.alloc(nb_element,
size_inverse_coords*size_inverse_coords,
type, ghost_type);
if(!interpolation_points_matrices.exists(type, ghost_type))
interpolation_points_matrices.alloc(nb_element,
nb_interpolation_points_per_elem * size_inverse_coords,
type, ghost_type);
Array<Real> & interp_inv_coord = interpolation_inverse_coordinates(type, ghost_type);
Array<Real> & interp_points_mat = interpolation_points_matrices(type, ghost_type);
Matrix<Real> quad_coord_matrix(size_inverse_coords, size_inverse_coords);
Array<Real>::const_iterator< Matrix<Real> > quad_coords_it =
quad_coordinates.begin_reinterpret(spatial_dimension,
nb_quad_per_element,
nb_element);
Array<Real>::const_iterator< Matrix<Real> > points_coords_begin =
interpolation_points_coordinates.begin_reinterpret(spatial_dimension,
nb_interpolation_points_per_elem,
interpolation_points_coordinates.getSize() / nb_interpolation_points_per_elem);
Array<Real>::iterator< Matrix<Real> > inv_quad_coord_it =
interp_inv_coord.begin(size_inverse_coords, size_inverse_coords);
Array<Real>::iterator< Matrix<Real> > inv_points_mat_it =
interp_points_mat.begin(nb_interpolation_points_per_elem, size_inverse_coords);
/// loop over the elements of the current material and element type
for (UInt el = 0; el < nb_element; ++el, ++inv_quad_coord_it,
++inv_points_mat_it, ++quad_coords_it) {
/// matrix containing the quadrature points coordinates
const Matrix<Real> & quad_coords = *quad_coords_it;
/// matrix to store the matrix inversion result
Matrix<Real> & inv_quad_coord_matrix = *inv_quad_coord_it;
/// insert the quad coordinates in a matrix compatible with the interpolation
buildElementalFieldInterpolationCoodinates<type>(quad_coords,
quad_coord_matrix);
/// invert the interpolation matrix
inv_quad_coord_matrix.inverse(quad_coord_matrix);
/// matrix containing the interpolation points coordinates
const Matrix<Real> & points_coords = points_coords_begin[elem_fil(el)];
/// matrix to store the interpolation points coordinates
/// compatible with these functions
Matrix<Real> & inv_points_coord_matrix = *inv_points_mat_it;
/// insert the quad coordinates in a matrix compatible with the interpolation
buildElementalFieldInterpolationCoodinates<type>(points_coords,
inv_points_coord_matrix);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void Material::interpolateStress(ByElementTypeReal & result,
const GhostType ghost_type) {
AKANTU_DEBUG_IN();
const Mesh & mesh = model->getFEM().getMesh();
Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type);
Mesh::type_iterator last = mesh.lastType(spatial_dimension, ghost_type);
for (; it != last; ++it) {
ElementType type = *it;
UInt nb_element = mesh.getNbElement(type, ghost_type);
if (nb_element == 0) continue;
Array<Real> & res = result(type, ghost_type);
#define INTERPOLATE_ELEMENTAL_FIELD(type) \
interpolateElementalField<type>(stress(type, ghost_type), \
res, \
ghost_type)
AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(INTERPOLATE_ELEMENTAL_FIELD);
#undef INTERPOLATE_ELEMENTAL_FIELD
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <ElementType type>
void Material::interpolateElementalField(const Array<Real> & field,
Array<Real> & result,
const GhostType ghost_type) {
AKANTU_DEBUG_IN();
Array<UInt> & elem_fil = element_filter(type, ghost_type);
UInt nb_element = elem_fil.getSize();
UInt nb_quad_per_element = model->getFEM().getNbQuadraturePoints(type, ghost_type);
UInt size_inverse_coords = getSizeElementalFieldInterpolationCoodinates<type>(ghost_type);
Matrix<Real> coefficients(nb_quad_per_element, field.getNbComponent());
const Array<Real> & interp_inv_coord = interpolation_inverse_coordinates(type,
ghost_type);
const Array<Real> & interp_points_coord = interpolation_points_matrices(type,
ghost_type);
UInt nb_interpolation_points_per_elem = interp_points_coord.getNbComponent() / size_inverse_coords;
Array<Real>::const_iterator< Matrix<Real> > field_it
= field.begin_reinterpret(field.getNbComponent(),
nb_quad_per_element,
nb_element);
Array<Real>::const_iterator< Matrix<Real> > interpolation_points_coordinates_it =
interp_points_coord.begin(nb_interpolation_points_per_elem, size_inverse_coords);
Array<Real>::iterator< Matrix<Real> > result_begin
= result.begin_reinterpret(field.getNbComponent(),
nb_interpolation_points_per_elem,
result.getSize() / nb_interpolation_points_per_elem);
Array<Real>::const_iterator< Matrix<Real> > inv_quad_coord_it =
interp_inv_coord.begin(size_inverse_coords, size_inverse_coords);
/// loop over the elements of the current material and element type
for (UInt el = 0; el < nb_element;
++el, ++field_it, ++inv_quad_coord_it, ++interpolation_points_coordinates_it) {
/**
* matrix containing the inversion of the quadrature points'
* coordinates
*/
const Matrix<Real> & inv_quad_coord_matrix = *inv_quad_coord_it;
/**
* multiply it by the field values over quadrature points to get
* the interpolation coefficients
*/
coefficients.mul<false, true>(inv_quad_coord_matrix, *field_it);
/// matrix containing the points' coordinates
const Matrix<Real> & coord = *interpolation_points_coordinates_it;
/// multiply the coordinates matrix by the coefficients matrix and store the result
result_begin[elem_fil(el)].mul<true, true>(coefficients, coord);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
const Array<Real> & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) const {
std::stringstream sstr;
std::string ghost_id = "";
if (ghost_type == _ghost) ghost_id = ":ghost";
sstr << getID() << ":" << vect_id << ":" << type << ghost_id;
ID fvect_id = sstr.str();
try {
return Memory::getArray<Real>(fvect_id);
} catch(debug::Exception & e) {
AKANTU_EXCEPTION("The material " << name << "(" <<getID() << ") does not contain a vector " << vect_id << "(" << fvect_id << ") [" << e << "]");
}
}
/* -------------------------------------------------------------------------- */
Array<Real> & Material::getArray(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) {
std::stringstream sstr;
std::string ghost_id = "";
if (ghost_type == _ghost) ghost_id = ":ghost";
sstr << getID() << ":" << vect_id << ":" << type << ghost_id;
ID fvect_id = sstr.str();
try {
return Memory::getArray<Real>(fvect_id);
} catch(debug::Exception & e) {
AKANTU_EXCEPTION("The material " << name << "(" << getID() << ") does not contain a vector " << vect_id << "(" << fvect_id << ") [" << e << "]");
}
}
/* -------------------------------------------------------------------------- */
const ByElementTypeArray<Real> & Material::getInternal(const ID & int_id) const {
std::map<ID, InternalField<Real> *>::const_iterator it = internal_vectors_real.find(getID() + ":" + int_id);
if(it == internal_vectors_real.end()) {
AKANTU_EXCEPTION("The material " << name << "(" << getID()
<< ") does not contain an internal "
<< int_id << " (" << (getID() + ":" + int_id) << ")");
}
return *it->second;
}
/* -------------------------------------------------------------------------- */
ByElementTypeArray<Real> & Material::getInternal(const ID & int_id) {
std::map<ID, InternalField<Real> *>::iterator it = internal_vectors_real.find(getID() + ":" + int_id);
if(it == internal_vectors_real.end()) {
AKANTU_EXCEPTION("The material " << name << "(" << getID()
<< ") does not contain an internal "
<< int_id << " (" << (getID() + ":" + int_id) << ")");
}
return *it->second;
}
/* -------------------------------------------------------------------------- */
void Material::printself(std::ostream & stream, int indent) const {
std::string space;
for(Int i = 0; i < indent; i++, space += AKANTU_INDENT);
std::string type = getID().substr(getID().find_last_of(":") + 1);
stream << space << "Material " << type << " [" << std::endl;
Parsable::printself(stream, indent);
stream << space << "]" << std::endl;
}
__END_AKANTU__
diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc
index 0f2be4d36..acfd27ab0 100644
--- a/src/model/solid_mechanics/solid_mechanics_model.cc
+++ b/src/model/solid_mechanics/solid_mechanics_model.cc
@@ -1,1881 +1,1891 @@
/**
* @file solid_mechanics_model.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Tue Jul 27 18:15:37 2010
*
* @brief Implementation of the SolidMechanicsModel class
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_math.hh"
#include "aka_common.hh"
#include "solid_mechanics_model.hh"
#include "integration_scheme_2nd_order.hh"
#include "element_group.hh"
#include "static_communicator.hh"
#include "dof_synchronizer.hh"
#include <cmath>
#ifdef AKANTU_USE_MUMPS
#include "solver_mumps.hh"
#endif
#ifdef AKANTU_USE_IOHELPER
# include "dumper_paraview.hh"
# include "dumper_iohelper_tmpl.hh"
# include "dumper_iohelper_tmpl_homogenizing_field.hh"
# include "dumper_iohelper_tmpl_material_internal_field.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), Dumpable(),
BoundaryCondition<SolidMechanicsModel>(),
time_step(NAN), f_m2a(1.0),
mass_matrix(NULL),
velocity_damping_matrix(NULL),
stiffness_matrix(NULL),
jacobian_matrix(NULL),
element_index_by_material("element index by material", id),
material_selector(new DefaultMaterialSelector(element_index_by_material)),
is_default_material_selector(true),
integrator(NULL),
increment_flag(false), solver(NULL),
synch_parallel(NULL),
are_materials_instantiated(false) {
AKANTU_DEBUG_IN();
createSynchronizerRegistry(this);
registerFEMObject<MyFEMType>("SolidMechanicsFEM", mesh, spatial_dimension);
this->displacement = NULL;
this->mass = NULL;
this->velocity = NULL;
this->acceleration = NULL;
this->force = NULL;
this->residual = NULL;
this->boundary = 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->registerDumper<DumperParaview>("paraview_all", id, true);
this->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;
if(solver) delete solver;
if(mass_matrix) delete mass_matrix;
if(velocity_damping_matrix) delete velocity_damping_matrix;
if(stiffness_matrix && stiffness_matrix != jacobian_matrix)
delete stiffness_matrix;
if(jacobian_matrix) delete jacobian_matrix;
if(dof_synchronizer) delete dof_synchronizer;
if(synch_parallel) delete synch_parallel;
if(is_default_material_selector) delete material_selector;
AKANTU_DEBUG_OUT();
}
void SolidMechanicsModel::setTimeStep(Real time_step) {
this->time_step = time_step;
#if defined(AKANTU_USE_IOHELPER)
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(std::string input_file, const ModelOptions & options) {
Model::initFull(input_file, 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(input_file != "") {
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::initFEMBoundary() {
FEM & fem_boundary = getFEMBoundary();
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 << ":boundary";
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));
boundary = &(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);
element_index_by_material.alloc(nb_element, 2, *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
getFEM().initShapeFunctions(_not_ghost);
getFEM().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)
(*boundary)((*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);
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->values;
Real * position_val = mesh.getNodes().values;
Real * displacement_val = displacement->values;
/// compute current_position = initial_position + displacement
memcpy(current_position_val, position_val, nb_nodes*spatial_dimension*sizeof(Real));
for (UInt n = 0; n < nb_nodes*spatial_dimension; ++n) {
*current_position_val++ += *displacement_val++;
}
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->values, force->values, nb_nodes*spatial_dimension*sizeof(Real));
// start synchronization
synch_registry->asynchronousSynchronize(_gst_smm_uv);
synch_registry->waitEndSynchronize(_gst_smm_uv);
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();
updateCurrentPosition();
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.savePreviousState(_not_ghost);
mat.savePreviousState(_ghost);
mat.computeAllStresses(_not_ghost);
if(mat.isFiniteDeformation())
mat.computeAllCauchyStresses(_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.savePreviousState(_not_ghost);
mat.savePreviousState(_ghost);
mat.computeAllStresses(_not_ghost);
if(mat.isFiniteDeformation())
mat.computeAllCauchyStresses(_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->values;
Real * accel_val = acceleration->values;
Real * res_val = residual->values;
bool * boundary_val = boundary->values;
for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) {
if(!(*boundary_val)) {
*res_val -= *accel_val * *mass_val /f_m2a;
}
boundary_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,
*boundary,
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> & boundary,
Real alpha) {
Real * A_val = A.storage();
Real * b_val = b.storage();
Real * x_val = x.storage();
bool * boundary_val = boundary.storage();
UInt nb_degrees_of_freedom = x.getSize() * x.getNbComponent();
for (UInt n = 0; n < nb_degrees_of_freedom; ++n) {
if(!(*boundary_val)) {
*x_val = alpha * (*b_val / *A_val);
}
x_val++;
A_val++;
b_val++;
boundary_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,
*boundary);
if(increment_flag) {
Real * inc_val = increment->values;
Real * dis_val = displacement->values;
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,
*boundary,
*increment_acceleration);
if(previous_displacement) previous_displacement->copy(*displacement);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::solveStep() {
AKANTU_DEBUG_IN();
this->explicitPred();
this->updateResidual();
this->updateAcceleration();
this->explicitCorr();
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) // 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";
jacobian_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension, _symmetric,
spatial_dimension, sstr.str(), memory_id);
jacobian_matrix->buildProfile(mesh, *dof_synchronizer);
if (!isExplicit()) {
delete stiffness_matrix;
std::stringstream sstr_sti; sstr_sti << id << ":stiffness_matrix";
stiffness_matrix = new SparseMatrix(*jacobian_matrix, sstr_sti.str(), memory_id);
}
#ifdef AKANTU_USE_MUMPS
std::stringstream sstr_solv; sstr_solv << id << ":solver";
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
}
/* -------------------------------------------------------------------------- */
/**
* 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(*boundary);
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();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::solveDynamic() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL,
"You should first initialize the implicit solver and assemble the stiffness matrix");
AKANTU_DEBUG_ASSERT(mass_matrix != NULL,
"You should first initialize the implicit solver and assemble the mass matrix");
// updateResidualInternal();
solve<NewmarkBeta::_displacement_corrector>(*increment);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::solveStatic() {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Solving Ku = f");
AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL,
"You should first initialize the implicit solver and assemble the stiffness matrix");
UInt nb_nodes = displacement->getSize();
UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes;
jacobian_matrix->copyContent(*stiffness_matrix);
jacobian_matrix->applyBoundary(*boundary);
increment->clear();
solver->setRHS(*residual);
solver->solve(*increment);
Real * increment_val = increment->storage();
Real * displacement_val = displacement->storage();
bool * boundary_val = boundary->storage();
for (UInt j = 0; j < nb_degree_of_freedom;
++j, ++displacement_val, ++increment_val, ++boundary_val) {
if (!(*boundary_val)) {
*displacement_val += *increment_val;
}
else {
*increment_val = 0.0;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::solveStatic(Array<bool> & boundary_normal, Array<Real> & EulerAngles) {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_INFO("Solving Ku = f");
AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL,
"You should first initialize the implicit solver and assemble the stiffness matrix");
UInt nb_nodes = displacement->getSize();
UInt nb_degree_of_freedom = displacement->getNbComponent();
// if(method != _static)
jacobian_matrix->copyContent(*stiffness_matrix);
Array<Real> * residual_rotated = new Array<Real > (nb_nodes, nb_degree_of_freedom, "residual_rotated");
//stiffness_matrix->saveMatrix("stiffness_original.out");
jacobian_matrix->applyBoundaryNormal(boundary_normal, EulerAngles, *residual, (*stiffness_matrix).getA(), *residual_rotated);
//jacobian_matrix->saveMatrix("stiffness_rotated_dir.out");
jacobian_matrix->applyBoundary(*boundary);
solver->setRHS(*residual_rotated);
delete residual_rotated;
if (!increment) setIncrementFlagOn();
solver->solve(*increment);
Matrix<Real> T(nb_degree_of_freedom, nb_degree_of_freedom);
Matrix<Real> small_rhs(nb_degree_of_freedom, nb_degree_of_freedom);
Matrix<Real> T_small_rhs(nb_degree_of_freedom, nb_degree_of_freedom);
Real * increment_val = increment->values;
Real * displacement_val = displacement->values;
bool * boundary_val = boundary->values;
for (UInt n = 0; n < nb_nodes; ++n) {
bool constrain_ij = false;
for (UInt j = 0; j < nb_degree_of_freedom; j++) {
if (boundary_normal(n, j)) {
constrain_ij = true;
break;
}
}
if (constrain_ij) {
if (nb_degree_of_freedom == 2) {
Real Theta = EulerAngles(n, 0);
T(0, 0) = cos(Theta);
T(0, 1) = -sin(Theta);
T(1, 1) = cos(Theta);
T(1, 0) = sin(Theta);
} else if (nb_degree_of_freedom == 3) {
Real Theta_x = EulerAngles(n, 0);
Real Theta_y = EulerAngles(n, 1);
Real Theta_z = EulerAngles(n, 2);
T(0, 0) = cos(Theta_y) * cos(Theta_z);
T(0, 1) = -cos(Theta_y) * sin(Theta_z);
T(0, 2) = sin(Theta_y);
T(1, 0) = cos(Theta_x) * sin(Theta_z) + cos(Theta_z) * sin(Theta_x) * sin(Theta_y);
T(1, 1) = cos(Theta_x) * cos(Theta_z) - sin(Theta_x) * sin(Theta_y) * sin(Theta_z);
T(1, 2) = -cos(Theta_y) * sin(Theta_x);
T(2, 0) = sin(Theta_x) * sin(Theta_z) - cos(Theta_x) * cos(Theta_z) * sin(Theta_y);
T(2, 1) = cos(Theta_z) * sin(Theta_x) + cos(Theta_x) * sin(Theta_y) * sin(Theta_z);
T(2, 2) = cos(Theta_x) * cos(Theta_y);
}
small_rhs.clear();
T_small_rhs.clear();
for (UInt j = 0; j < nb_degree_of_freedom; j++)
if(!(boundary_normal(n,j)) )
small_rhs(j,j)=increment_val[j];
T_small_rhs.mul<true, false>(T,small_rhs);
for (UInt j = 0; j < nb_degree_of_freedom; j++){
if(!(boundary_normal(n,j))){
for (UInt k = 0; k < nb_degree_of_freedom; k++)
displacement_val[k]+=T_small_rhs(k,j);
}
}
displacement_val += nb_degree_of_freedom;
boundary_val += nb_degree_of_freedom;
increment_val += nb_degree_of_freedom;
} else {
for (UInt j = 0; j < nb_degree_of_freedom; j++, ++displacement_val, ++increment_val, ++boundary_val) {
if (!(*boundary_val)) {
*displacement_val += *increment_val;
}
}
}
}
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 * boundary_val = boundary->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(!(*boundary_val) && is_local_node) {
norm[0] += *increment_val * *increment_val;
norm[1] += *displacement_val * *displacement_val;
}
boundary_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");
AKANTU_DEBUG_OUT();
error = norm[0] / norm[1];
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->values;
bool * boundary_val = boundary->values;
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(!(*boundary_val)) {
norm += *residual_val * *residual_val;
}
boundary_val++;
residual_val++;
}
} else {
boundary_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);
}
/* -------------------------------------------------------------------------- */
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,
*boundary);
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::implicitCorr() {
AKANTU_DEBUG_IN();
if(method == _implicit_dynamic) {
integrator->integrationSchemeCorrDispl(time_step,
*displacement,
*velocity,
*acceleration,
*boundary,
*increment);
} else {
UInt nb_nodes = displacement->getSize();
UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes;
Real * incr_val = increment->values;
Real * disp_val = displacement->values;
bool * boun_val = boundary->values;
for (UInt j = 0; j < nb_degree_of_freedom; ++j, ++disp_val, ++incr_val, ++boun_val){
*disp_val += (1. - *boun_val) * *incr_val;
*incr_val *= (1. - *boun_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->values;
Real * disp_val = displacement->values;
Real * prev_disp_val = previous_displacement->values;
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>::iterator< Vector<UInt> > eibm =
element_index_by_material(*it, ghost_type).begin(2);
Array<Real> X(0, nb_nodes_per_element*spatial_dimension);
FEM::extractNodalToElementField(mesh, *current_position,
X, *it, _not_ghost);
Array<Real>::iterator< Matrix<Real> > X_el =
X.begin(spatial_dimension, nb_nodes_per_element);
for (UInt el = 0; el < nb_element; ++el, ++X_el, ++eibm) {
elem.element = el;
Real el_size = getFEM().getElementInradius(*X_el, *it);
Real el_dt = mat_val[(*eibm)(0)]->getStableTimeStep(el_size, elem);
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->values;
Real * mass_val = mass->values;
for (UInt n = 0; n < nb_nodes; ++n) {
Real mv2 = 0;
bool is_local_node = mesh.isLocalOrMasterNode(n);
bool is_not_pbc_slave_node = !getIsPBCSlaveNode(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 = getFEM().getNbQuadraturePoints(type);
Array<Real> vel_on_quad(nb_quadrature_points, spatial_dimension);
Array<UInt> filter_element(1, 1, index);
getFEM().interpolateOnQuadraturePoints(*velocity, vel_on_quad,
spatial_dimension,
type, _not_ghost,
filter_element);
Array<Real>::iterator< Vector<Real> > vit = vel_on_quad.begin(spatial_dimension);
Array<Real>::iterator< Vector<Real> > vend = vel_on_quad.end(spatial_dimension);
Vector<Real> rho_v2(nb_quadrature_points);
Real rho = materials[element_index_by_material(type)(index, 0)]->getRho();
for (UInt q = 0; vit != vend; ++vit, ++q) {
rho_v2(q) = rho * vit->dot(*vit);
}
AKANTU_DEBUG_OUT();
return .5*getFEM().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 = boundary->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 = !getIsPBCSlaveNode(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;
Vector<UInt> mat = element_index_by_material(type, _not_ghost).begin(2)[index];
Real energy = materials[mat(0)]->getEnergy(energy_id, type, mat(1));
AKANTU_DEBUG_OUT();
return energy;
}
/* -------------------------------------------------------------------------- */
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(boundary ) boundary ->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) {
delete stiffness_matrix;
delete jacobian_matrix;
delete solver;
SolverOptions solver_options;
initImplicit((method == _implicit_dynamic), solver_options);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onElementsAdded(const Array<Element> & element_list,
const NewElementsEvent & event) {
AKANTU_DEBUG_IN();
getFEM().initShapeFunctions(_not_ghost);
getFEM().initShapeFunctions(_ghost);
Array<Element>::const_iterator<Element> it = element_list.begin();
Array<Element>::const_iterator<Element> end = element_list.end();
/// \todo have rules to choose the correct material
UInt mat_id = 0;
UInt * mat_id_vect = NULL;
try {
const NewMaterialElementsEvent & event_mat = dynamic_cast<const NewMaterialElementsEvent &>(event);
mat_id_vect = event_mat.getMaterialList().storage();
} catch(...) { }
for (UInt el = 0; it != end; ++it, ++el) {
const Element & elem = *it;
if(mat_id_vect) mat_id = mat_id_vect[el];
else mat_id = (*material_selector)(elem);
Material & mat = *materials[mat_id];
UInt mat_index = mat.addElement(elem.type, elem.element, elem.ghost_type);
Vector<UInt> id(2);
id[0] = mat_id; id[1] = mat_index;
if(!element_index_by_material.exists(elem.type, elem.ghost_type))
element_index_by_material.alloc(0, 2, elem.type, elem.ghost_type);
element_index_by_material(elem.type, elem.ghost_type).push_back(id);
}
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) AKANTU_DEBUG_TO_IMPLEMENT();
assembleMassLumped();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::onElementsRemoved(__attribute__((unused)) const Array<Element> & element_list,
const ByElementTypeUInt & new_numbering,
const RemovedElementsEvent & event) {
// MeshUtils::purifyMesh(mesh);
getFEM().initShapeFunctions(_not_ghost);
getFEM().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::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(boundary ) mesh.removeNodesFromArray(*boundary , 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();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::addDumpFieldToDumper(const std::string & dumper_name,
const std::string & field_id) {
#ifdef AKANTU_USE_IOHELPER
#define ADD_FIELD(field, type) \
internalAddDumpFieldToDumper(dumper_name, \
BOOST_PP_STRINGIZE(field), \
new DumperIOHelper::NodalField<type>(*field))
if(field_id == "displacement") { ADD_FIELD(displacement, Real); }
else if(field_id == "mass" ) { ADD_FIELD(mass , Real); }
else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); }
else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); }
else if(field_id == "force" ) { ADD_FIELD(force , Real); }
else if(field_id == "residual" ) { ADD_FIELD(residual , Real); }
else if(field_id == "boundary" ) { ADD_FIELD(boundary , bool); }
else if(field_id == "increment" ) { ADD_FIELD(increment , Real); }
else if(field_id == "partitions" ) {
internalAddDumpFieldToDumper(dumper_name,
field_id,
new DumperIOHelper::ElementPartitionField<>(mesh,
spatial_dimension,
_not_ghost,
_ek_regular));
}
else if(field_id == "element_index_by_material") {
internalAddDumpFieldToDumper(dumper_name,
field_id,
new DumperIOHelper::ElementalField<UInt>(element_index_by_material,
spatial_dimension,
_not_ghost,
_ek_regular));
} else {
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_id);
}
if (is_internal) {
internalAddDumpFieldToDumper
(dumper_name,
field_id,
new DumperIOHelper::HomogenizedField<Real,
DumperIOHelper::InternalMaterialField>
(*this,
field_id,
spatial_dimension,
_not_ghost,
_ek_regular));
} else {
try {
internalAddDumpFieldToDumper
(dumper_name,
field_id,
new DumperIOHelper::ElementalField<UInt>(mesh.getData<UInt>(field_id),
spatial_dimension,
_not_ghost,
_ek_regular));
} catch (...) {}
+ try {
+ internalAddDumpFieldToDumper
+ (dumper_name,
+ field_id,
+ new DumperIOHelper::ElementalField<Real>(mesh.getData<Real>(field_id),
+ spatial_dimension,
+ _not_ghost,
+ _ek_regular));
+ } catch (...) {}
}
}
#undef ADD_FIELD
#endif
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::addDumpGroupField(const std::string & field_id,
const std::string & group_name) {
ElementGroup & group = mesh.getElementGroup(group_name);
this->addDumpGroupFieldToDumper(group.getDefaultDumperName(),
field_id,
group_name);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::addDumpGroupFieldToDumper(const std::string & dumper_name,
const std::string & field_id,
const std::string & group_name) {
#ifdef AKANTU_USE_IOHELPER
ElementGroup & group = mesh.getElementGroup(group_name);
UInt dim = group.getDimension();
const Array<UInt> * nodal_filter = &(group.getNodes());
const ByElementTypeArray<UInt> * elemental_filter = &(group.getElements());
#define ADD_FIELD(field, type) \
group.addDumpFieldExternalToDumper(dumper_name, \
BOOST_PP_STRINGIZE(field), \
new DumperIOHelper::NodalField<type, true>(*field, 0, 0, nodal_filter))
if(field_id == "displacement") { ADD_FIELD(displacement, Real); }
else if(field_id == "mass" ) { ADD_FIELD(mass , Real); }
else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); }
else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); }
else if(field_id == "force" ) { ADD_FIELD(force , Real); }
else if(field_id == "residual" ) { ADD_FIELD(residual , Real); }
else if(field_id == "boundary" ) { ADD_FIELD(boundary , bool); }
else if(field_id == "increment" ) { ADD_FIELD(increment , Real); }
else if(field_id == "partitions" ) {
group.addDumpFieldExternalToDumper(dumper_name,
field_id,
new DumperIOHelper::ElementPartitionField<true>(mesh,
dim,
_not_ghost,
_ek_regular,
elemental_filter));
} else if(field_id == "element_index_by_material") {
group.addDumpFieldExternalToDumper(dumper_name,
field_id,
new DumperIOHelper::ElementalField<UInt, Vector, true>(element_index_by_material,
dim,
_not_ghost,
_ek_regular,
elemental_filter));
} else {
typedef DumperIOHelper::HomogenizedField<Real,
DumperIOHelper::InternalMaterialField,
DumperIOHelper::internal_material_field_iterator,
DumperIOHelper::AvgHomogenizingFunctor,
Vector,
true> Field;
group.addDumpFieldExternalToDumper(dumper_name,
field_id,
new Field(*this,
field_id,
dim,
_not_ghost,
_ek_regular,
elemental_filter));
}
#undef ADD_FIELD
#endif
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::removeDumpGroupField(const std::string & field_id,
const std::string & group_name) {
ElementGroup & group = mesh.getElementGroup(group_name);
this->removeDumpGroupFieldFromDumper(group.getDefaultDumperName(),
field_id,
group_name);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::removeDumpGroupFieldFromDumper(const std::string & dumper_name,
const std::string & field_id,
const std::string & group_name) {
ElementGroup & group = mesh.getElementGroup(group_name);
group.removeDumpFieldFromDumper(dumper_name, field_id);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::addDumpFieldVectorToDumper(const std::string & dumper_name,
const std::string & field_id) {
#ifdef AKANTU_USE_IOHELPER
#define ADD_FIELD(field, type) \
DumperIOHelper::Field * f = \
new DumperIOHelper::NodalField<type>(*field); \
f->setPadding(3); \
internalAddDumpFieldToDumper(dumper_name, BOOST_PP_STRINGIZE(field), f)
if(field_id == "displacement") { ADD_FIELD(displacement, Real); }
else if(field_id == "mass" ) { ADD_FIELD(mass , Real); }
else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); }
else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); }
else if(field_id == "force" ) { ADD_FIELD(force , Real); }
else if(field_id == "residual" ) { ADD_FIELD(residual , Real); }
else {
typedef DumperIOHelper::HomogenizedField<Real,
DumperIOHelper::InternalMaterialField> Field;
Field * field = new Field(*this,
field_id,
spatial_dimension,
spatial_dimension,
_not_ghost,
_ek_regular);
field->setPadding(3);
internalAddDumpFieldToDumper(dumper_name, field_id, field);
}
#undef ADD_FIELD
#endif
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::addDumpGroupFieldVector(const std::string & field_id,
const std::string & group_name) {
ElementGroup & group = mesh.getElementGroup(group_name);
this->addDumpGroupFieldVectorToDumper(group.getDefaultDumperName(),
field_id,
group_name);
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::addDumpGroupFieldVectorToDumper(const std::string & dumper_name,
const std::string & field_id,
const std::string & group_name) {
#ifdef AKANTU_USE_IOHELPER
ElementGroup & group = mesh.getElementGroup(group_name);
UInt dim = group.getDimension();
const Array<UInt> * nodal_filter = &(group.getNodes());
const ByElementTypeArray<UInt> * elemental_filter = &(group.getElements());
#define ADD_FIELD(field, type) \
DumperIOHelper::Field * f = \
new DumperIOHelper::NodalField<type, true>(*field, 0, 0, nodal_filter); \
f->setPadding(3); \
group.addDumpFieldExternalToDumper(dumper_name, BOOST_PP_STRINGIZE(field), f)
if(field_id == "displacement") { ADD_FIELD(displacement, Real); }
else if(field_id == "mass" ) { ADD_FIELD(mass , Real); }
else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); }
else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); }
else if(field_id == "force" ) { ADD_FIELD(force , Real); }
else if(field_id == "residual" ) { ADD_FIELD(residual , Real); }
else {
typedef DumperIOHelper::HomogenizedField<Real,
DumperIOHelper::InternalMaterialField,
DumperIOHelper::internal_material_field_iterator,
DumperIOHelper::AvgHomogenizingFunctor,
Vector,
true> Field;
Field * field = new Field(*this,
field_id,
spatial_dimension,
dim,
_not_ghost,
_ek_regular,
elemental_filter);
field->setPadding(3);
group.addDumpFieldExternalToDumper(dumper_name, field_id, field);
}
#undef ADD_FIELD
#endif
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModel::addDumpFieldTensorToDumper(const std::string & dumper_name,
const std::string & field_id) {
#ifdef AKANTU_USE_IOHELPER
if(field_id == "stress") {
typedef DumperIOHelper::HomogenizedField<Real,
DumperIOHelper::InternalMaterialField,
DumperIOHelper::material_stress_field_iterator,
DumperIOHelper::AvgHomogenizingFunctor,
Matrix> Field;
Field * field = new Field(*this,
field_id,
spatial_dimension,
spatial_dimension,
_not_ghost,
_ek_regular);
field->setPadding(3, 3);
internalAddDumpFieldToDumper(dumper_name, field_id, field);
} else if(field_id == "strain") {
typedef DumperIOHelper::HomogenizedField<Real,
DumperIOHelper::InternalMaterialField,
DumperIOHelper::material_strain_field_iterator,
DumperIOHelper::AvgHomogenizingFunctor,
Matrix> Field;
Field * field = new Field(*this,
field_id,
spatial_dimension,
spatial_dimension,
_not_ghost,
_ek_regular);
field->setPadding(3, 3);
internalAddDumpFieldToDumper(dumper_name, field_id, field);
} else {
typedef DumperIOHelper::HomogenizedField<Real,
DumperIOHelper::InternalMaterialField,
DumperIOHelper::internal_material_field_iterator,
DumperIOHelper::AvgHomogenizingFunctor,
Matrix> Field;
Field * field = new Field(*this,
field_id,
spatial_dimension,
spatial_dimension,
_not_ghost,
_ek_regular);
field->setPadding(3, 3);
internalAddDumpFieldToDumper(dumper_name, field_id, field);
}
#endif
}
/* -------------------------------------------------------------------------- */
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;
getFEM().printself(stream, indent + 2);
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << " + nodals information [" << std::endl;
displacement->printself(stream, indent + 2);
mass ->printself(stream, indent + 2);
velocity ->printself(stream, indent + 2);
acceleration->printself(stream, indent + 2);
force ->printself(stream, indent + 2);
residual ->printself(stream, indent + 2);
boundary ->printself(stream, indent + 2);
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << " + connectivity type information [" << std::endl;
element_index_by_material.printself(stream, indent + 2);
stream << space << AKANTU_INDENT << "]" << std::endl;
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
__END_AKANTU__
diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc
index 1507e389b..66333c339 100644
--- a/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc
+++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive.cc
@@ -1,1001 +1,807 @@
/**
* @file solid_mechanics_model_cohesive.cc
*
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date Tue May 08 13:01:18 2012
*
* @brief Solid mechanics model for cohesive elements
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include <algorithm>
#include "shape_cohesive.hh"
#include "solid_mechanics_model_cohesive.hh"
#include "material_cohesive.hh"
#ifdef AKANTU_USE_IOHELPER
# include "dumper_paraview.hh"
# include "dumper_iohelper_tmpl.hh"
# include "dumper_iohelper_tmpl_homogenizing_field.hh"
# include "dumper_iohelper_tmpl_material_internal_field.hh"
#endif
/* -------------------------------------------------------------------------- */
__BEGIN_AKANTU__
const SolidMechanicsModelCohesiveOptions default_solid_mechanics_model_cohesive_options(_explicit_lumped_mass,
false,
false,
true,
true);
/* -------------------------------------------------------------------------- */
SolidMechanicsModelCohesive::SolidMechanicsModelCohesive(Mesh & mesh,
UInt dim,
const ID & id,
const MemoryID & memory_id) :
SolidMechanicsModel(mesh, dim, id, memory_id),
mesh_facets(mesh.initMeshFacets("mesh_facets")),
tangents("tangents", id),
stress_on_facet("stress_on_facet", id),
facet_stress("facet_stress", id),
- fragment_mass(0, spatial_dimension, "fragment_mass"),
- fragment_velocity(0, spatial_dimension, "fragment_velocity"),
- fragment_center(0, spatial_dimension, "fragment_center"),
facet_material("facet_material", id),
inserter(mesh, mesh_facets) {
AKANTU_DEBUG_IN();
facet_generated = false;
#if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT)
facet_synchronizer = NULL;
facet_stress_synchronizer = NULL;
cohesive_distributed_synchronizer = NULL;
global_connectivity = NULL;
#endif
delete material_selector;
material_selector = new DefaultMaterialCohesiveSelector(*this);
#if defined(AKANTU_USE_IOHELPER)
this->registerDumper<DumperParaview>("cohesive_elements", id);
this->addDumpMeshToDumper("cohesive_elements",
mesh, spatial_dimension, _not_ghost, _ek_cohesive);
#endif
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
SolidMechanicsModelCohesive::~SolidMechanicsModelCohesive() {
AKANTU_DEBUG_IN();
#if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT)
delete cohesive_distributed_synchronizer;
delete facet_synchronizer;
delete facet_stress_synchronizer;
#endif
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::setTimeStep(Real time_step) {
SolidMechanicsModel::setTimeStep(time_step);
#if defined(AKANTU_USE_IOHELPER)
getDumper("cohesive_elements").setTimeStep(time_step);
#endif
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::initFull(std::string material_file,
const ModelOptions & options) {
AKANTU_DEBUG_IN();
const SolidMechanicsModelCohesiveOptions & smmc_options =
dynamic_cast<const SolidMechanicsModelCohesiveOptions &>(options);
this->stress_interpolation = smmc_options.stress_interpolation;
this->is_extrinsic = smmc_options.extrinsic;
if (is_extrinsic) {
if (!facet_generated)
MeshUtils::buildAllFacets(mesh, mesh_facets);
}
SolidMechanicsModel::initFull(material_file, options);
if (is_extrinsic)
initAutomaticInsertion();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::initMaterials() {
AKANTU_DEBUG_IN();
/// init element inserter
inserter.init();
#if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT)
if (facet_synchronizer != NULL)
inserter.initParallel(facet_synchronizer);
if (facet_stress_synchronizer != NULL) {
DataAccessor * data_accessor = this;
const ByElementTypeUInt & rank_to_element = synch_parallel->getPrankToElement();
facet_stress_synchronizer->updateFacetStressSynchronizer(inserter,
rank_to_element,
*data_accessor);
}
#endif
// make sure the material are instantiated
if(!are_materials_instantiated) instantiateMaterials();
/// find the first cohesive material
UInt cohesive_index = 0;
while ((dynamic_cast<MaterialCohesive *>(materials[cohesive_index]) == NULL)
&& cohesive_index <= materials.size())
++cohesive_index;
AKANTU_DEBUG_ASSERT(cohesive_index != materials.size(),
"No cohesive materials in the material input file");
material_selector->setFallback(cohesive_index);
// set the facet information in the material in case of dynamic insertion
mesh_facets.initByElementTypeArray(facet_material, 1, spatial_dimension - 1);
if (is_extrinsic) {
Element element;
for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) {
element.ghost_type = *gt;
Mesh::type_iterator first = mesh_facets.firstType(spatial_dimension - 1, *gt);
Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1, *gt);
for(;first != last; ++first) {
element.type = *first;
Array<UInt> & f_material = facet_material(*first, *gt);
UInt nb_element = mesh_facets.getNbElement(*first, *gt);
f_material.resize(nb_element);
f_material.set(cohesive_index);
for (UInt el = 0; el < nb_element; ++el) {
element.element = el;
UInt mat_index = (*material_selector)(element);
f_material(el) = mat_index;
MaterialCohesive & mat = dynamic_cast<MaterialCohesive &>(*materials[mat_index]);
mat.addFacet(element);
}
}
}
} else {
for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) {
Mesh::type_iterator first = mesh.firstType(spatial_dimension, *gt, _ek_cohesive);
Mesh::type_iterator last = mesh.lastType(spatial_dimension, *gt, _ek_cohesive);
for(;first != last; ++first) {
Array<UInt> & el_id_by_mat = element_index_by_material(*first, *gt);
Vector<UInt> el_mat(2); el_mat(0) = cohesive_index; el_mat(1) = 0;
el_id_by_mat.set(el_mat);
}
}
}
SolidMechanicsModel::initMaterials();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
/**
* Initialize the model,basically it pre-compute the shapes, shapes derivatives
* and jacobian
*
*/
void SolidMechanicsModelCohesive::initModel() {
AKANTU_DEBUG_IN();
SolidMechanicsModel::initModel();
registerFEMObject<MyFEMCohesiveType>("CohesiveFEM", mesh, spatial_dimension);
/// add cohesive type connectivity
ElementType type = _not_defined;
for (ghost_type_t::iterator gt = ghost_type_t::begin();
gt != ghost_type_t::end(); ++gt) {
GhostType type_ghost = *gt;
Mesh::type_iterator it = mesh.firstType(spatial_dimension, type_ghost);
Mesh::type_iterator last = mesh.lastType(spatial_dimension, type_ghost);
for (; it != last; ++it) {
const Array<UInt> & connectivity = mesh.getConnectivity(*it, type_ghost);
if (connectivity.getSize() != 0) {
type = *it;
ElementType type_facet = Mesh::getFacetType(type);
ElementType type_cohesive = FEM::getCohesiveElementType(type_facet);
mesh.addConnectivityType(type_cohesive, type_ghost);
}
}
}
AKANTU_DEBUG_ASSERT(type != _not_defined, "No elements in the mesh");
getFEM("CohesiveFEM").initShapeFunctions(_not_ghost);
getFEM("CohesiveFEM").initShapeFunctions(_ghost);
registerFEMObject<MyFEMType>("FacetsFEM", mesh_facets, spatial_dimension - 1);
if (is_extrinsic) {
getFEM("FacetsFEM").initShapeFunctions(_not_ghost);
getFEM("FacetsFEM").initShapeFunctions(_ghost);
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::initAutomaticInsertion() {
AKANTU_DEBUG_IN();
mesh_facets.initByElementTypeArray(facet_stress,
2 * spatial_dimension * spatial_dimension,
spatial_dimension - 1);
resizeFacetStress();
/// compute normals on facets
computeNormals();
/// allocate stress_on_facet to store element stress on facets
mesh.initByElementTypeArray(stress_on_facet,
spatial_dimension * spatial_dimension,
spatial_dimension);
Mesh::type_iterator it = mesh.firstType(spatial_dimension);
Mesh::type_iterator last = mesh.lastType(spatial_dimension);
for (; it != last; ++it) {
ElementType type = *it;
UInt nb_element = mesh.getNbElement(type);
UInt nb_facet_per_elem = Mesh::getNbFacetsPerElement(type);
ElementType type_facet = Mesh::getFacetType(type);
UInt nb_quad_per_facet = getFEM("FacetsFEM").getNbQuadraturePoints(type_facet);
stress_on_facet(type).resize(nb_quad_per_facet * nb_facet_per_elem * nb_element);
}
if (stress_interpolation)
initStressInterpolation();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::updateAutomaticInsertion() {
AKANTU_DEBUG_IN();
inserter.limitCheckFacets();
#if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT)
if (facet_stress_synchronizer != NULL) {
DataAccessor * data_accessor = this;
const ByElementTypeUInt & rank_to_element = synch_parallel->getPrankToElement();
facet_stress_synchronizer->updateFacetStressSynchronizer(inserter,
rank_to_element,
*data_accessor);
}
#endif
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::initStressInterpolation() {
/// compute quadrature points coordinates on facets
Array<Real> & position = mesh.getNodes();
ByElementTypeReal quad_facets("quad_facets", id);
mesh_facets.initByElementTypeArray(quad_facets,
spatial_dimension, spatial_dimension - 1);
- for (ghost_type_t::iterator gt = ghost_type_t::begin();
- gt != ghost_type_t::end(); ++gt) {
-
- GhostType facet_gt = *gt;
- Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, facet_gt);
- Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1, facet_gt);
-
- for (; it != last; ++it) {
- ElementType facet_type = *it;
-
- UInt nb_quad_per_facet =
- getFEM("FacetsFEM").getNbQuadraturePoints(facet_type);
-
- UInt nb_facet = mesh_facets.getNbElement(facet_type, facet_gt);
- UInt nb_tot_quad = nb_quad_per_facet * nb_facet;
-
- Array<Real> & quad_f = quad_facets(facet_type, facet_gt);
- quad_f.resize(nb_tot_quad);
-
- getFEM("FacetsFEM").interpolateOnQuadraturePoints(position,
- quad_f,
- spatial_dimension,
- facet_type,
- facet_gt);
- }
- }
+ getFEM("FacetsFEM").interpolateOnQuadraturePoints(position, quad_facets);
/// compute elements quadrature point positions and build
/// element-facet quadrature points data structure
ByElementTypeReal elements_quad_facets("elements_quad_facets", id);
mesh.initByElementTypeArray(elements_quad_facets,
spatial_dimension,
spatial_dimension);
for (ghost_type_t::iterator gt = ghost_type_t::begin();
gt != ghost_type_t::end(); ++gt) {
GhostType elem_gt = *gt;
Mesh::type_iterator it = mesh.firstType(spatial_dimension, elem_gt);
Mesh::type_iterator last = mesh.lastType(spatial_dimension, elem_gt);
for (; it != last; ++it) {
ElementType type = *it;
UInt nb_element = mesh.getNbElement(type, elem_gt);
if (nb_element == 0) continue;
/// compute elements' quadrature points and list of facet
/// quadrature points positions by element
Array<Element> & facet_to_element = mesh_facets.getSubelementToElement(type,
elem_gt);
UInt nb_facet_per_elem = facet_to_element.getNbComponent();
Array<Real> & el_q_facet = elements_quad_facets(type, elem_gt);
ElementType facet_type = Mesh::getFacetType(type);
UInt nb_quad_per_facet =
getFEM("FacetsFEM").getNbQuadraturePoints(facet_type);
el_q_facet.resize(nb_element * nb_facet_per_elem * nb_quad_per_facet);
for (UInt el = 0; el < nb_element; ++el) {
for (UInt f = 0; f < nb_facet_per_elem; ++f) {
Element global_facet_elem = facet_to_element(el, f);
UInt global_facet = global_facet_elem.element;
GhostType facet_gt = global_facet_elem.ghost_type;
const Array<Real> & quad_f = quad_facets(facet_type, facet_gt);
for (UInt q = 0; q < nb_quad_per_facet; ++q) {
for (UInt s = 0; s < spatial_dimension; ++s) {
el_q_facet(el * nb_facet_per_elem * nb_quad_per_facet
+ f * nb_quad_per_facet + q, s)
= quad_f(global_facet * nb_quad_per_facet + q, s);
}
}
}
}
}
}
/// loop over non cohesive materials
for (UInt m = 0; m < materials.size(); ++m) {
try {
MaterialCohesive & mat __attribute__((unused)) =
dynamic_cast<MaterialCohesive &>(*materials[m]);
} catch(std::bad_cast&) {
/// initialize the interpolation function
materials[m]->initElementalFieldInterpolation(elements_quad_facets);
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::updateResidual(bool need_initialize) {
AKANTU_DEBUG_IN();
if (need_initialize) initializeUpdateResidualData();
// f -= fint
std::vector<Material *>::iterator mat_it;
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
try {
MaterialCohesive & mat = dynamic_cast<MaterialCohesive &>(**mat_it);
mat.computeTraction(_not_ghost);
} catch (std::bad_cast & bce) { }
}
SolidMechanicsModel::updateResidual(false);
for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) {
try {
MaterialCohesive & mat = dynamic_cast<MaterialCohesive &>(**mat_it);
mat.computeEnergies();
} catch (std::bad_cast & bce) { }
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::computeNormals() {
AKANTU_DEBUG_IN();
getFEM("FacetsFEM").computeNormalsOnControlPoints(_not_ghost);
/**
* @todo store tangents while computing normals instead of
* recomputing them as follows:
*/
/* ------------------------------------------------------------------------ */
UInt tangent_components = spatial_dimension * (spatial_dimension - 1);
mesh_facets.initByElementTypeArray(tangents,
tangent_components,
spatial_dimension - 1);
Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1);
Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1);
for (; it != last; ++it) {
ElementType facet_type = *it;
const Array<Real> & normals =
getFEM("FacetsFEM").getNormalsOnQuadPoints(facet_type);
UInt nb_quad = normals.getSize();
Array<Real> & tang = tangents(facet_type);
tang.resize(nb_quad);
Real * normal_it = normals.storage();
Real * tangent_it = tang.storage();
/// compute first tangent
for (UInt q = 0; q < nb_quad; ++q) {
/// if normal is orthogonal to xy plane, arbitrarly define tangent
if ( Math::are_float_equal(Math::norm2(normal_it), 0) )
tangent_it[0] = 1;
else
Math::normal2(normal_it, tangent_it);
normal_it += spatial_dimension;
tangent_it += tangent_components;
}
/// compute second tangent (3D case)
if (spatial_dimension == 3) {
normal_it = normals.storage();
tangent_it = tang.storage();
for (UInt q = 0; q < nb_quad; ++q) {
Math::normal3(normal_it, tangent_it, tangent_it + spatial_dimension);
normal_it += spatial_dimension;
tangent_it += tangent_components;
}
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::averageStressOnFacets(UInt material_index) {
AKANTU_DEBUG_IN();
Mesh::type_iterator it = mesh.firstType(spatial_dimension);
Mesh::type_iterator last = mesh.lastType(spatial_dimension);
/// loop over element type
for (; it != last; ++it) {
ElementType type = *it;
UInt nb_element = mesh.getNbElement(type);
UInt nb_quad_per_element = getFEM().getNbQuadraturePoints(type);
const Array<Real> & stress = materials[material_index]->getStress(type);
Array<Real> & s_on_facet = stress_on_facet(type);
UInt nb_facet_per_elem = Mesh::getNbFacetsPerElement(type);
ElementType type_facet = Mesh::getFacetType(type);
UInt nb_quad_per_facet = getFEM("FacetsFEM").getNbQuadraturePoints(type_facet);
UInt nb_facet_quad_per_elem = nb_quad_per_facet * nb_facet_per_elem;
Array<Real>::const_iterator<Matrix<Real> > stress_it
= stress.begin(spatial_dimension, spatial_dimension);
Array<Real>::iterator<Matrix<Real> > s_on_facet_it
= s_on_facet.begin(spatial_dimension, spatial_dimension);
Matrix<Real> average_stress(spatial_dimension, spatial_dimension);
for (UInt el = 0; el < nb_element; ++el) {
/// compute average stress
average_stress.clear();
for (UInt q = 0; q < nb_quad_per_element; ++q, ++stress_it)
average_stress += *stress_it;
average_stress /= nb_quad_per_element;
/// store average stress
for (UInt q = 0; q < nb_facet_quad_per_elem; ++q, ++s_on_facet_it)
*s_on_facet_it = average_stress;
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::fillStressOnFacet() {
AKANTU_DEBUG_IN();
UInt sp2 = spatial_dimension * spatial_dimension;
UInt sp4 = sp2 * 2;
/// loop over materials
for (UInt m = 0; m < materials.size(); ++m) {
try {
MaterialCohesive & mat __attribute__((unused)) =
dynamic_cast<MaterialCohesive &>(*materials[m]);
} catch(std::bad_cast&) {
if (stress_interpolation)
/// interpolate stress on facet quadrature points positions
materials[m]->interpolateStress(stress_on_facet);
else
averageStressOnFacets(m);
GhostType ghost_type = _not_ghost;
Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type);
Mesh::type_iterator last = mesh.lastType(spatial_dimension, ghost_type);
/// loop over element type
for (; it != last; ++it) {
ElementType type = *it;
UInt nb_element = mesh.getNbElement(type, ghost_type);
if (nb_element == 0) continue;
Array<Real> & stress_on_f = stress_on_facet(type, ghost_type);
/// store the interpolated stresses on the facet_stress vector
Array<Element> & facet_to_element =
mesh_facets.getSubelementToElement(type, ghost_type);
UInt nb_facet_per_elem = facet_to_element.getNbComponent();
Array<Element>::iterator<Vector<Element> > facet_to_el_it =
facet_to_element.begin(nb_facet_per_elem);
Array<Real>::iterator< Matrix<Real> > stress_on_f_it =
stress_on_f.begin(spatial_dimension, spatial_dimension);
ElementType facet_type = _not_defined;
GhostType facet_gt = _casper;
Array<std::vector<Element> > * element_to_facet = NULL;
Array<Real> * f_stress = NULL;
Array<bool> * f_check = NULL;
UInt nb_quad_per_facet = 0;
UInt element_rank = 0;
Element current_el(type, 0, ghost_type);
for (UInt el = 0; el < nb_element; ++el, ++facet_to_el_it) {
current_el.element = el;
for (UInt f = 0; f < nb_facet_per_elem; ++f) {
Element global_facet_elem = (*facet_to_el_it)(f);
UInt global_facet = global_facet_elem.element;
if (facet_type != global_facet_elem.type ||
facet_gt != global_facet_elem.ghost_type) {
facet_type = global_facet_elem.type;
facet_gt = global_facet_elem.ghost_type;
element_to_facet =
&( mesh_facets.getElementToSubelement(facet_type, facet_gt) );
f_stress = &( facet_stress(facet_type, facet_gt) );
nb_quad_per_facet =
getFEM("FacetsFEM").getNbQuadraturePoints(facet_type, facet_gt);
f_check = &( inserter.getCheckFacets(facet_type, facet_gt) );
}
if (!(*f_check)(global_facet)) {
stress_on_f_it += nb_quad_per_facet;
continue;
}
for (UInt q = 0; q < nb_quad_per_facet; ++q, ++stress_on_f_it) {
element_rank = (*element_to_facet)(global_facet)[0] != current_el;
Matrix<Real> facet_stress_local(f_stress->storage()
+ (global_facet * nb_quad_per_facet + q) * sp4
+ element_rank * sp2,
spatial_dimension,
spatial_dimension);
facet_stress_local = *stress_on_f_it;
}
}
}
}
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::checkCohesiveStress() {
AKANTU_DEBUG_IN();
fillStressOnFacet();
#if defined(AKANTU_DEBUG_TOOLS)
debug::element_manager.printData(debug::_dm_model_cohesive,
"Interpolated stresses before",
facet_stress);
#endif
synch_registry->synchronize(_gst_smmc_facets_stress);
#if defined(AKANTU_DEBUG_TOOLS)
debug::element_manager.printData(debug::_dm_model_cohesive,
"Interpolated stresses",
facet_stress);
#endif
for (UInt m = 0; m < materials.size(); ++m) {
try {
MaterialCohesive & mat_cohesive = dynamic_cast<MaterialCohesive &>(*materials[m]);
/// check which not ghost cohesive elements are to be created
mat_cohesive.checkInsertion();
} catch(std::bad_cast&) { }
}
/// communicate data among processors
synch_registry->synchronize(_gst_smmc_facets);
/// insert cohesive elements
inserter.insertExtrinsicElements();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::onElementsAdded(const Array<Element> & element_list,
const NewElementsEvent & event) {
AKANTU_DEBUG_IN();
SolidMechanicsModel::onElementsAdded(element_list, event);
/// update shape functions
getFEM("CohesiveFEM").initShapeFunctions(_not_ghost);
getFEM("CohesiveFEM").initShapeFunctions(_ghost);
if(is_extrinsic) {
getFEM("FacetsFEM").initShapeFunctions(_not_ghost);
getFEM("FacetsFEM").initShapeFunctions(_ghost);
}
#if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT)
updateFacetSynchronizers();
#endif
resizeFacetStress();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::onNodesAdded(const Array<UInt> & doubled_nodes,
__attribute__((unused)) const NewNodesEvent & event) {
AKANTU_DEBUG_IN();
UInt nb_new_nodes = doubled_nodes.getSize();
Array<UInt> nodes_list(nb_new_nodes);
for (UInt n = 0; n < nb_new_nodes; ++n)
nodes_list(n) = doubled_nodes(n, 1);
SolidMechanicsModel::onNodesAdded(nodes_list, event);
for (UInt n = 0; n < nb_new_nodes; ++n) {
UInt old_node = doubled_nodes(n, 0);
UInt new_node = doubled_nodes(n, 1);
for (UInt dim = 0; dim < spatial_dimension; ++dim) {
(*displacement)(new_node, dim) = (*displacement)(old_node, dim);
(*velocity) (new_node, dim) = (*velocity) (old_node, dim);
(*acceleration)(new_node, dim) = (*acceleration)(old_node, dim);
(*boundary) (new_node, dim) = (*boundary) (old_node, dim);
if (current_position)
(*current_position)(new_node, dim) = (*current_position)(old_node, dim);
if (increment_acceleration)
(*increment_acceleration)(new_node, dim)
= (*increment_acceleration)(old_node, dim);
if (increment)
(*increment)(new_node, dim) = (*increment)(old_node, dim);
- }
- }
-
- AKANTU_DEBUG_OUT();
-}
-
-/* -------------------------------------------------------------------------- */
-class CohesiveElementFilter : public GroupManager::ClusteringFilter {
-public:
- CohesiveElementFilter(const SolidMechanicsModelCohesive & model,
- const Real max_damage = 1.) :
- model(model), is_unbroken(max_damage) {}
-
- bool operator() (const Element & el) const {
- if (el.kind == _ek_regular)
- return true;
-
- const Array<UInt> & el_id_by_mat = model.getElementIndexByMaterial(el.type,
- el.ghost_type);
-
- const MaterialCohesive & mat
- = static_cast<const MaterialCohesive &>
- (model.getMaterial(el_id_by_mat(el.element, 0)));
-
- UInt el_index = el_id_by_mat(el.element, 1);
- UInt nb_quad_per_element
- = model.getFEM("CohesiveFEM").getNbQuadraturePoints(el.type, el.ghost_type);
-
- Vector<Real> damage
- = mat.getDamage(el.type, el.ghost_type).begin(nb_quad_per_element)[el_index];
-
- UInt unbroken_quads = std::count_if(damage.storage(),
- damage.storage() + nb_quad_per_element,
- is_unbroken);
-
- if (unbroken_quads > 0)
- return true;
-
- return false;
- }
-
-private:
-
- struct IsUnbrokenFunctor {
- IsUnbrokenFunctor(const Real & max_damage) : max_damage(max_damage) {}
- bool operator() (const Real & x) {return x < max_damage;}
- const Real max_damage;
- };
-
- const SolidMechanicsModelCohesive & model;
- const IsUnbrokenFunctor is_unbroken;
-};
-
-/* -------------------------------------------------------------------------- */
-void SolidMechanicsModelCohesive::buildFragments() {
- AKANTU_DEBUG_IN();
-
-#if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT)
- if (cohesive_distributed_synchronizer) {
- cohesive_distributed_synchronizer->computeBufferSize(*this, _gst_smmc_damage);
- cohesive_distributed_synchronizer->asynchronousSynchronize(*this, _gst_smmc_damage);
- cohesive_distributed_synchronizer->waitEndSynchronize(*this, _gst_smmc_damage);
- }
-#endif
- nb_fragment = mesh.createClusters(spatial_dimension, "fragment",
- CohesiveElementFilter(*this),
- synch_parallel,
- &mesh_facets);
-
- AKANTU_DEBUG_OUT();
-}
-
-/* -------------------------------------------------------------------------- */
-void SolidMechanicsModelCohesive::computeFragmentsMassVelocity() {
- AKANTU_DEBUG_IN();
-
- fragment_mass.resize(nb_fragment);
- fragment_mass.clear();
-
- fragment_velocity.resize(nb_fragment);
- fragment_velocity.clear();
-
- fragment_center.resize(nb_fragment);
- fragment_center.clear();
-
- Array<Real>::const_iterator< Vector<Real> > it_mass =
- mass->begin(spatial_dimension);
- Array<Real>::const_iterator< Vector<Real> > it_velocity =
- velocity->begin(spatial_dimension);
- Array<Real>::const_iterator< Vector<Real> > it_position =
- mesh.getNodes().begin(spatial_dimension);
-
- Array<Real>::iterator< Vector<Real> > it_frag_mass =
- fragment_mass.begin(spatial_dimension);
- Array<Real>::iterator< Vector<Real> > it_frag_velocity =
- fragment_velocity.begin(spatial_dimension);
- Array<Real>::iterator< Vector<Real> > it_frag_center =
- fragment_center.begin(spatial_dimension);
-
-
- /// compute data for local fragments
- for (UInt frag = 0; frag < nb_fragment; ++frag) {
- std::stringstream sstr;
- sstr << "fragment_" << frag;
-
- GroupManager::const_element_group_iterator eg_it
- = mesh.element_group_find(sstr.str());
-
- if (eg_it == mesh.element_group_end()) continue;
-
- Vector<Real> & frag_mass = it_frag_mass[frag];
- Vector<Real> & frag_velocity = it_frag_velocity[frag];
- Vector<Real> & frag_center = it_frag_center[frag];
-
- const NodeGroup & node_group = eg_it->second->getNodeGroup();
- NodeGroup::const_node_iterator it = node_group.begin();
- NodeGroup::const_node_iterator end = node_group.end();
-
- for (; it != end; ++it) {
- UInt node = *it;
- if (mesh.isLocalOrMasterNode(node)) {
- const Vector<Real> & node_mass = it_mass[node];
- const Vector<Real> & node_velocity = it_velocity[node];
- const Vector<Real> & node_position = it_position[node];
-
- frag_mass += node_mass;
-
- for (UInt s = 0; s < spatial_dimension; ++s) {
- frag_velocity(s)
- += node_mass(s) * node_velocity(s);
-
- frag_center(s)
- += node_mass(s) * node_position(s);
- }
- }
- }
- }
-
-#if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT)
- /// reduce all data
- StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator();
- comm.allReduce(fragment_mass.storage() , nb_fragment * spatial_dimension, _so_sum);
- comm.allReduce(fragment_velocity.storage(), nb_fragment * spatial_dimension, _so_sum);
- comm.allReduce(fragment_center.storage() , nb_fragment * spatial_dimension, _so_sum);
-#endif
-
- /// compute averages
- it_frag_mass = fragment_mass.begin(spatial_dimension);
- it_frag_velocity = fragment_velocity.begin(spatial_dimension);
- it_frag_center = fragment_center.begin(spatial_dimension);
-
- for (UInt frag = 0;
- frag < nb_fragment;
- ++frag, ++it_frag_mass, ++it_frag_velocity, ++it_frag_center) {
-
- for (UInt s = 0; s < spatial_dimension; ++s) {
- (*it_frag_velocity)(s) /= (*it_frag_mass)(s);
- (*it_frag_center)(s) /= (*it_frag_mass)(s);
+ if (previous_displacement)
+ (*previous_displacement)(new_node, dim)
+ = (*previous_displacement)(old_node, dim);
}
}
AKANTU_DEBUG_OUT();
}
-/* -------------------------------------------------------------------------- */
-void SolidMechanicsModelCohesive::computeFragmentsData() {
- AKANTU_DEBUG_IN();
-
- buildFragments();
- computeFragmentsMassVelocity();
-
- AKANTU_DEBUG_OUT();
-}
-
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::printself(std::ostream & stream, int indent) const {
std::string space;
for(Int i = 0; i < indent; i++, space += AKANTU_INDENT);
stream << space << "SolidMechanicsModelCohesive [" << std::endl;
SolidMechanicsModel::printself(stream, indent + 1);
stream << space << "]" << std::endl;
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::resizeFacetStress() {
AKANTU_DEBUG_IN();
for (ghost_type_t::iterator gt = ghost_type_t::begin();
gt != ghost_type_t::end();
++gt) {
GhostType ghost_type = *gt;
Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, ghost_type);
Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1, ghost_type);
for(; it != end; ++it) {
ElementType type = *it;
UInt nb_facet = mesh_facets.getNbElement(type, ghost_type);
UInt nb_quadrature_points =
getFEM("FacetsFEM").getNbQuadraturePoints(type, ghost_type);
UInt new_size = nb_facet * nb_quadrature_points;
facet_stress(type, ghost_type).resize(new_size);
}
}
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
void SolidMechanicsModelCohesive::addDumpFieldToDumper(const std::string & dumper_name,
const std::string & field_id) {
AKANTU_DEBUG_IN();
if (dumper_name == "cohesive_elements") {
if (field_id == "damage") {
internalAddDumpFieldToDumper
("cohesive_elements",
field_id, new DumperIOHelper::
HomogenizedField<Real,
DumperIOHelper::InternalMaterialField>(*this,
field_id,
spatial_dimension,
_not_ghost,
_ek_cohesive));
}
} else {
SolidMechanicsModel::addDumpFieldToDumper(dumper_name, field_id);
}
AKANTU_DEBUG_OUT();
}
-
__END_AKANTU__
diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive.hh
index c76a56d90..a815c4d6e 100644
--- a/src/model/solid_mechanics/solid_mechanics_model_cohesive.hh
+++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive.hh
@@ -1,315 +1,282 @@
/**
* @file solid_mechanics_model_cohesive.hh
*
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date Tue May 08 13:01:18 2012
*
* @brief Solid mechanics model for cohesive elements
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
+#ifndef __AKANTU_SOLID_MECHANICS_MODEL_COHESIVE_HH__
+#define __AKANTU_SOLID_MECHANICS_MODEL_COHESIVE_HH__
+
#include "solid_mechanics_model.hh"
#include "cohesive_element_inserter.hh"
#if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT)
# include "facet_synchronizer.hh"
# include "facet_stress_synchronizer.hh"
#endif
/* -------------------------------------------------------------------------- */
-#ifndef __AKANTU_SOLID_MECHANICS_MODEL_COHESIVE_HH__
-#define __AKANTU_SOLID_MECHANICS_MODEL_COHESIVE_HH__
-
__BEGIN_AKANTU__
/* -------------------------------------------------------------------------- */
struct SolidMechanicsModelCohesiveOptions : public SolidMechanicsModelOptions {
SolidMechanicsModelCohesiveOptions(AnalysisMethod analysis_method = _explicit_lumped_mass,
bool extrinsic = false,
bool no_init_materials = false,
bool init_facet_filter = true,
bool stress_interpolation = true) :
SolidMechanicsModelOptions(analysis_method, no_init_materials),
extrinsic(extrinsic), init_facet_filter(init_facet_filter),
stress_interpolation(stress_interpolation) {}
bool extrinsic;
bool init_facet_filter;
bool stress_interpolation;
};
extern const SolidMechanicsModelCohesiveOptions default_solid_mechanics_model_cohesive_options;
/* -------------------------------------------------------------------------- */
/* Solid Mechanics Model for Cohesive elements */
/* -------------------------------------------------------------------------- */
class SolidMechanicsModelCohesive : public SolidMechanicsModel {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
class NewCohesiveNodesEvent : public NewNodesEvent {
public:
AKANTU_GET_MACRO_NOT_CONST(OldNodesList, old_nodes, Array<UInt> &);
AKANTU_GET_MACRO(OldNodesList, old_nodes, const Array<UInt> &);
protected:
Array<UInt> old_nodes;
};
typedef FEMTemplate<IntegratorGauss, ShapeLagrange, _ek_cohesive> MyFEMCohesiveType;
SolidMechanicsModelCohesive(Mesh & mesh,
UInt spatial_dimension = _all_dimensions,
const ID & id = "solid_mechanics_model_cohesive",
const MemoryID & memory_id = 0);
virtual ~SolidMechanicsModelCohesive();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// set the value of the time step
void setTimeStep(Real time_step);
/// assemble the residual for the explicit scheme
virtual void updateResidual(bool need_initialize = true);
/// function to print the contain of the class
virtual void printself(std::ostream & stream, int indent = 0) const;
/// function to perform a stress check on each facet and insert
/// cohesive elements if needed
void checkCohesiveStress();
/// initialize the cohesive model
void initFull(std::string material_file,
const ModelOptions & options = default_solid_mechanics_model_cohesive_options);
/// initialize the model
void initModel();
/// initialize cohesive material
void initMaterials();
- /// build fragments
- void buildFragments();
-
- /// compute fragments' mass and velocity
- void computeFragmentsMassVelocity();
-
- /// build fragments and compute their data (mass, velocity..)
- void computeFragmentsData();
-
/// init facet filters for cohesive materials
void initFacetFilter();
/// limit the cohesive element insertion to a given area
void enableFacetsCheckOnArea(const Array<Real> & limits);
/// update automatic insertion after a change in the element inserter
void updateAutomaticInsertion();
private:
/// initialize completely the model for extrinsic elements
void initAutomaticInsertion();
/// initialize stress interpolation
void initStressInterpolation();
/// compute facets' normals
void computeNormals();
/// resize facet stress
void resizeFacetStress();
/// init facets_check array
void initFacetsCheck();
/// fill stress_on_facet
void fillStressOnFacet();
/// compute average stress on elements
void averageStressOnFacets(UInt material_index);
/* ------------------------------------------------------------------------ */
/* Mesh Event Handler inherited members */
/* ------------------------------------------------------------------------ */
protected:
virtual void onNodesAdded (const Array<UInt> & nodes_list,
const NewNodesEvent & event);
virtual void onElementsAdded (const Array<Element> & nodes_list,
const NewElementsEvent & event);
/* ------------------------------------------------------------------------ */
/* Dumpable interface */
/* ------------------------------------------------------------------------ */
public:
virtual void addDumpFieldToDumper(const std::string & dumper_name,
const std::string & field_id);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/// get facet mesh
AKANTU_GET_MACRO(MeshFacets, mesh_facets, const Mesh &);
/// get stress on facets vector
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(StressOnFacets, facet_stress, Real);
- /// get number of fragments
- AKANTU_GET_MACRO(NbFragment, nb_fragment, UInt);
-
- /// get mass for each fragment
- AKANTU_GET_MACRO(FragmentsMass, fragment_mass, const Array<Real> &);
-
- /// get average velocity for each fragment
- AKANTU_GET_MACRO(FragmentsVelocity, fragment_velocity, const Array<Real> &);
-
- /// get center of mass coordinates for each fragment
- AKANTU_GET_MACRO(FragmentsCenter, fragment_center, const Array<Real> &);
-
/// get facet material
AKANTU_GET_MACRO_BY_ELEMENT_TYPE(FacetMaterial, facet_material, UInt);
/// get facet material
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(FacetMaterial, facet_material, UInt);
/// get facet material
AKANTU_GET_MACRO(FacetMaterial, facet_material, const ByElementTypeArray<UInt> &);
- /// THIS HAS TO BE CHANGED
+ /// @todo THIS HAS TO BE CHANGED
AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Tangents, tangents, Real);
/// get element inserter
AKANTU_GET_MACRO_NOT_CONST(ElementInserter, inserter, CohesiveElementInserter &);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
/// mesh containing facets and their data structures
Mesh & mesh_facets;
/// @todo store tangents when normals are computed:
ByElementTypeReal tangents;
/// list of stresses on facet quadrature points for every element
ByElementTypeReal stress_on_facet;
/// stress on facets on the two sides by quadrature point
ByElementTypeReal facet_stress;
- /// number of fragments
- UInt nb_fragment;
-
- /// mass for each fragment
- Array<Real> fragment_mass;
-
- /// average velocity for each fragment
- Array<Real> fragment_velocity;
-
- /// center of mass coordinates for each element
- Array<Real> fragment_center;
-
/// flag to know if facets have been generated
bool facet_generated;
/// material to use if a cohesive element is created on a facet
ByElementTypeUInt facet_material;
/// stress interpolation flag
bool stress_interpolation;
bool is_extrinsic;
/// cohesive element inserter
CohesiveElementInserter inserter;
#if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT)
#include "solid_mechanics_model_cohesive_parallel.hh"
#endif
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
#if defined (AKANTU_PARALLEL_COHESIVE_ELEMENT)
# include "solid_mechanics_model_cohesive_inline_impl.cc"
#endif
/* -------------------------------------------------------------------------- */
class DefaultMaterialCohesiveSelector : public DefaultMaterialSelector {
public:
DefaultMaterialCohesiveSelector(const SolidMechanicsModelCohesive & model) :
DefaultMaterialSelector(model.getElementIndexByMaterial()),
facet_material(model.getFacetMaterial()),
mesh(model.getMesh()),
mesh_facets(model.getMeshFacets()) { }
inline virtual UInt operator()(const Element & element) {
if(Mesh::getKind(element.type) == _ek_cohesive) {
try {
const Array<Element> & cohesive_el_to_facet
= mesh_facets.getSubelementToElement(element.type, element.ghost_type);
bool third_dimension = (mesh.getSpatialDimension() == 3);
const Element & facet = cohesive_el_to_facet(element.element, third_dimension);
return facet_material(facet.type, facet.ghost_type)(facet.element);
} catch(...) {
return MaterialSelector::operator()(element);
}
} else if (Mesh::getSpatialDimension(element.type) == mesh.getSpatialDimension() - 1) {
return facet_material(element.type, element.ghost_type)(element.element);
} else {
return DefaultMaterialSelector::operator()(element);
}
}
private:
const ByElementTypeUInt & facet_material;
const Mesh & mesh;
const Mesh & mesh_facets;
};
/// standard output stream operator
inline std::ostream & operator <<(std::ostream & stream, const SolidMechanicsModelCohesive & _this)
{
_this.printself(stream);
return stream;
}
__END_AKANTU__
#endif /* __AKANTU_SOLID_MECHANICS_MODEL_COHESIVE_HH__ */
diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive_parallel.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive_parallel.hh
index ee9d4074e..68d47fbf3 100644
--- a/src/model/solid_mechanics/solid_mechanics_model_cohesive_parallel.hh
+++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive_parallel.hh
@@ -1,104 +1,114 @@
/**
* @file solid_mechanics_model_cohesive_parallel.hh
*
* @author Marco Vocialta <marco.vocialta@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date Fri Jun 14 17:09:09 2013
*
* @brief Include of the class members for parallelism
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
/// register the tags associated with the parallel synchronizer for
/// cohesive elements
void initParallel(MeshPartition * partition,
DataAccessor * data_accessor = NULL,
bool extrinsic = false);
protected:
void fillSynchronizeNormals(ByElementTypeReal & facet_normals,
GhostType ghost_type);
void synchronizeGhostFacets();
void updateFacetSynchronizers();
void fillGlobalConnectivity(ByElementTypeUInt & global_connectivity,
GhostType ghost_type,
bool just_init);
+/* ------------------------------------------------------------------------ */
+/* Accessors */
+/* ------------------------------------------------------------------------ */
+public:
+
+/// get cohesive elements synchronizer
+AKANTU_GET_MACRO(CohesiveSynchronizer,
+ cohesive_distributed_synchronizer,
+ const DistributedSynchronizer *);
+
/* ------------------------------------------------------------------------ */
/* Data Accessor inherited members */
/* ------------------------------------------------------------------------ */
public:
inline UInt getNbQuadsForFacetCheck(const Array<Element> & elements) const;
inline UInt getNbNodesPerElementList(const Array<Element> & elements) const;
inline virtual UInt getNbDataForElements(const Array<Element> & elements,
SynchronizationTag tag) const;
inline virtual void packElementData(CommunicationBuffer & buffer,
const Array<Element> & elements,
SynchronizationTag tag) const;
inline virtual void unpackElementData(CommunicationBuffer & buffer,
const Array<Element> & elements,
SynchronizationTag tag);
template<typename T>
inline void packFacetStressDataHelper(const ByElementTypeArray<T> & data_to_pack,
CommunicationBuffer & buffer,
const Array<Element> & elements) const;
template<typename T>
inline void unpackFacetStressDataHelper(ByElementTypeArray<T> & data_to_unpack,
CommunicationBuffer & buffer,
const Array<Element> & elements) const;
template<typename T, bool pack_helper>
inline void packUnpackFacetStressDataHelper(ByElementTypeArray<T> & data_to_pack,
CommunicationBuffer & buffer,
const Array<Element> & element) const;
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
/// facet synchronizer
FacetSynchronizer * facet_synchronizer;
/// facet stress synchronizer
FacetStressSynchronizer * facet_stress_synchronizer;
/// cohesive elements synchronizer
DistributedSynchronizer * cohesive_distributed_synchronizer;
/// global connectivity
ByElementTypeUInt * global_connectivity;
diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_buildfragments/test_cohesive_buildfragments.cc b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_buildfragments/test_cohesive_buildfragments.cc
index 26c4fa724..85f84a0e0 100644
--- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_buildfragments/test_cohesive_buildfragments.cc
+++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_buildfragments/test_cohesive_buildfragments.cc
@@ -1,202 +1,205 @@
/**
* @file test_cohesive_buildfragments.cc
*
* @author Marco Vocialta <marco.vocialta@epfl.ch>
*
* @date Wed Jun 13 11:29:49 2012
*
* @brief Test for cohesive elements
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include <limits>
#include <fstream>
#include <iostream>
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "mesh.hh"
#include "mesh_utils.hh"
#include "solid_mechanics_model_cohesive.hh"
#include "material_cohesive.hh"
+#include "fragment_manager.hh"
+
//#include "io_helper.hh"
/* -------------------------------------------------------------------------- */
using namespace akantu;
int main(int argc, char *argv[]) {
initialize(argc, argv);
debug::setDebugLevel(dblWarning);
const UInt spatial_dimension = 2;
const UInt max_steps = 200;
Real strain_rate = 1.e5;
ElementType type = _quadrangle_4;
Real L = 0.03;
Real theoretical_mass = L * L/20. * 2500;
ElementType type_facet = Mesh::getFacetType(type);
ElementType type_cohesive = FEM::getCohesiveElementType(type_facet);
Mesh mesh(spatial_dimension);
mesh.read("mesh.msh");
mesh.createGroupsFromMeshData<std::string>("physical_names");
SolidMechanicsModelCohesive model(mesh);
/// model initialization
model.initFull("material.dat",
SolidMechanicsModelCohesiveOptions(_explicit_lumped_mass, true));
Real time_step = model.getStableTimeStep()*0.05;
model.setTimeStep(time_step);
// std::cout << "Time step: " << time_step << std::endl;
Real disp_increment = strain_rate * L / 2. * time_step;
model.assembleMassLumped();
Array<Real> & velocity = model.getVelocity();
const Array<Real> & position = mesh.getNodes();
UInt nb_nodes = mesh.getNbNodes();
/// initial conditions
for (UInt n = 0; n < nb_nodes; ++n)
velocity(n, 0) = strain_rate * position(n, 0);
/// boundary conditions
model.applyBC(BC::Dirichlet::FixedValue(0, BC::_x), "Left_side");
model.applyBC(BC::Dirichlet::FixedValue(0, BC::_x), "Right_side");
model.updateResidual();
// model.setBaseName("extrinsic");
// model.addDumpFieldVector("displacement");
// model.addDumpField("velocity" );
// model.addDumpField("acceleration");
// model.addDumpField("residual" );
// model.addDumpField("stress");
// model.addDumpField("strain");
// model.dump();
UInt cohesive_index = 1;
UInt nb_quad_per_facet = model.getFEM("FacetsFEM").getNbQuadraturePoints(type_facet);
MaterialCohesive & mat_cohesive
= dynamic_cast<MaterialCohesive&>(model.getMaterial(cohesive_index));
const Array<Real> & damage = mat_cohesive.getDamage(type_cohesive);
- const Array<Real> & fragment_mass = model.getFragmentsMass();
+ FragmentManager fragment_manager(model, false);
+ const Array<Real> & fragment_mass = fragment_manager.getMass();
/// Main loop
for (UInt s = 1; s <= max_steps; ++s) {
model.checkCohesiveStress();
model.explicitPred();
model.updateResidual();
model.updateAcceleration();
model.explicitCorr();
/// apply boundary conditions
model.applyBC(BC::Dirichlet::IncrementValue(-disp_increment, BC::_x), "Left_side");
model.applyBC(BC::Dirichlet::IncrementValue( disp_increment, BC::_x), "Right_side");
if(s % 1 == 0) {
// model.dump();
std::cout << "passing step " << s << "/" << max_steps << std::endl;
- model.computeFragmentsData();
+ fragment_manager.computeAllData();
/// check number of fragments
- UInt nb_fragment_num = model.getNbFragment();
+ UInt nb_fragment_num = fragment_manager.getNbFragment();
UInt nb_cohesive_elements = mesh.getNbElement(type_cohesive);
UInt nb_fragment = 1;
for (UInt el = 0; el < nb_cohesive_elements; ++el) {
UInt q = 0;
while (q < nb_quad_per_facet &&
Math::are_float_equal(damage(el*nb_quad_per_facet + q), 1)) ++q;
if (q == nb_quad_per_facet) {
++nb_fragment;
}
}
if (nb_fragment != nb_fragment_num) {
std::cout << "The number of fragments is wrong!" << std::endl;
return EXIT_FAILURE;
}
/// check mass computation
Real total_mass = 0.;
for (UInt frag = 0; frag < nb_fragment_num; ++frag) {
total_mass += fragment_mass(frag);
}
if (!Math::are_float_equal(theoretical_mass, total_mass)) {
std::cout << "The fragments' mass is wrong!" << std::endl;
return EXIT_FAILURE;
}
}
}
/// check velocities
- UInt nb_fragment = model.getNbFragment();
- const Array<Real> & fragment_velocity = model.getFragmentsVelocity();
- const Array<Real> & fragment_center = model.getFragmentsCenter();
+ UInt nb_fragment = fragment_manager.getNbFragment();
+ const Array<Real> & fragment_velocity = fragment_manager.getVelocity();
+ const Array<Real> & fragment_center = fragment_manager.getCenterOfMass();
Real fragment_length = L / nb_fragment;
Real initial_position = -L / 2. + fragment_length / 2.;
for (UInt frag = 0; frag < nb_fragment; ++frag) {
Real theoretical_center = initial_position + fragment_length * frag;
if (!Math::are_float_equal(fragment_center(frag, 0), theoretical_center)) {
std::cout << "The fragments' center is wrong!" << std::endl;
return EXIT_FAILURE;
}
Real initial_vel = fragment_center(frag, 0) * strain_rate;
Math::setTolerance(100);
if (!Math::are_float_equal(fragment_velocity(frag), initial_vel)) {
std::cout << "The fragments' velocity is wrong!" << std::endl;
return EXIT_FAILURE;
}
}
finalize();
std::cout << "OK: test_cohesive_buildfragments was passed!" << std::endl;
return EXIT_SUCCESS;
}
diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/material.dat b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/material.dat
index b2e9636dc..e2f11daf6 100644
--- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/material.dat
+++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_extrinsic/material.dat
@@ -1,22 +1,23 @@
-material elastic [
+material plastic [
name = steel
rho = 1e3 # density
E = 1e3 # young's modulus
nu = 0.001 # poisson's ratio
+ finite_deformation = 1
]
material cohesive_linear [
name = cohesive
beta = 1
G_cI = 100
G_cII = 100
sigma_c = 100
]
material cohesive_linear [
name = cohesive2
beta = 1
G_cI = 100
G_cII = 100
sigma_c = 1000
]
diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_parallel_cohesive/test_cohesive_parallel_buildfragments/test_cohesive_parallel_buildfragments.cc b/test/test_model/test_solid_mechanics_model/test_cohesive/test_parallel_cohesive/test_cohesive_parallel_buildfragments/test_cohesive_parallel_buildfragments.cc
index 4ec4affd5..2c9ea53e6 100644
--- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_parallel_cohesive/test_cohesive_parallel_buildfragments/test_cohesive_parallel_buildfragments.cc
+++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_parallel_cohesive/test_cohesive_parallel_buildfragments/test_cohesive_parallel_buildfragments.cc
@@ -1,323 +1,327 @@
/**
* @file test_cohesive_parallel_buildfragments.cc
* @author Marco Vocialta <marco.vocialta@epfl.ch>
* @date Tue Dec 17 12:03:03 2013
*
* @brief Test to build fragments in parallel
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#include <limits>
#include <fstream>
#include <iostream>
/* -------------------------------------------------------------------------- */
#include "aka_common.hh"
#include "mesh.hh"
#include "mesh_utils.hh"
#include "solid_mechanics_model_cohesive.hh"
#include "material_cohesive.hh"
+#include "fragment_manager.hh"
-#ifdef AKANTU_USE_IOHELPER
-# include "dumper_paraview.hh"
-# include "dumper_iohelper_tmpl.hh"
-# include "dumper_iohelper_tmpl_homogenizing_field.hh"
-# include "dumper_iohelper_tmpl_material_internal_field.hh"
-#endif
//#include "io_helper.hh"
/* -------------------------------------------------------------------------- */
using namespace akantu;
void verticalInsertionLimit(SolidMechanicsModelCohesive &);
void displaceElements(SolidMechanicsModelCohesive &, const Real, const Real);
int main(int argc, char *argv[]) {
initialize(argc, argv);
debug::setDebugLevel(dblWarning);
const UInt spatial_dimension = 3;
const UInt total_nb_fragment = 50;
Mesh mesh(spatial_dimension);
StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator();
Int psize = comm.getNbProc();
Int prank = comm.whoAmI();
akantu::MeshPartition * partition = NULL;
if(prank == 0) {
// Read the mesh
mesh.read("mesh.msh");
/// partition the mesh
MeshUtils::purifyMesh(mesh);
partition = new MeshPartitionScotch(mesh, spatial_dimension);
partition->partitionate(psize);
}
SolidMechanicsModelCohesive model(mesh);
model.initParallel(partition, NULL, true);
delete partition;
/// model initialization
model.initFull("material.dat",
SolidMechanicsModelCohesiveOptions(_explicit_lumped_mass, true));
mesh.computeBoundingBox();
Real L = mesh.getXMax() - mesh.getXMin();
Real h = mesh.getYMax() - mesh.getYMin();
Real rho = model.getMaterial("bulk").getParam<Real>("rho");
Real theoretical_mass = L * h * h * rho;
Real frag_theo_mass = theoretical_mass / total_nb_fragment;
- model.computeFragmentsData();
- mesh.createMeshDataFromClusters("fragment");
+ FragmentManager fragment_manager(model);
+
+ fragment_manager.computeAllData();
model.setBaseName("extrinsic");
model.addDumpFieldVector("displacement");
model.addDumpField("velocity" );
model.addDumpField("stress");
model.addDumpField("partitions");
- model.addDumpField("fragment");
+ model.addDumpField("fragments");
+ model.addDumpField("fragments mass");
+ model.addDumpField("moments of inertia");
model.dump();
model.setBaseNameToDumper("cohesive_elements", "cohesive_elements_test");
model.addDumpFieldVectorToDumper("cohesive_elements", "displacement");
model.addDumpFieldToDumper("cohesive_elements", "damage");
model.dump("cohesive_elements");
/// set check facets
verticalInsertionLimit(model);
model.assembleMassLumped();
model.synchronizeBoundaries();
/// impose initial displacement
Array<Real> & displacement = model.getDisplacement();
Array<Real> & velocity = model.getVelocity();
const Array<Real> & position = mesh.getNodes();
UInt nb_nodes = mesh.getNbNodes();
for (UInt n = 0; n < nb_nodes; ++n) {
displacement(n, 0) = position(n, 0) * 0.1;
velocity(n, 0) = position(n, 0);
}
model.updateResidual();
model.checkCohesiveStress();
model.dump();
model.dump("cohesive_elements");
- const Array<Real> & fragment_mass = model.getFragmentsMass();
+ const Array<Real> & fragment_mass = fragment_manager.getMass();
+ const Array<Real> & fragment_center = fragment_manager.getCenterOfMass();
Real el_size = L / total_nb_fragment;
Real lim = -L/2 + el_size * 0.99;
/// displace one fragment each time
for (UInt frag = 1; frag <= total_nb_fragment; ++frag) {
if (prank == 0)
std::cout << "Generating fragment: " << frag << std::endl;
- model.computeFragmentsData();
- mesh.createMeshDataFromClusters("fragment");
+ fragment_manager.computeAllData();
/// check number of fragments
- UInt nb_fragment_num = model.getNbFragment();
+ UInt nb_fragment_num = fragment_manager.getNbFragment();
if (nb_fragment_num != frag) {
AKANTU_DEBUG_ERROR("The number of fragments is wrong! Numerical: " << nb_fragment_num << " Theoretical: " << frag);
return EXIT_FAILURE;
}
/// check mass computation
if (frag < total_nb_fragment) {
Real total_mass = 0.;
UInt small_fragments = 0;
for (UInt f = 0; f < nb_fragment_num; ++f) {
if (Math::are_float_equal(fragment_mass(f, 0), frag_theo_mass)) {
+
+ /// check center of mass
+ if (fragment_center(f, 0) > (L * frag / total_nb_fragment - L / 2)) {
+ std::cout << "Fragment center is wrong!" << std::endl;
+ return EXIT_FAILURE;
+ }
+
++small_fragments;
total_mass += frag_theo_mass;
}
}
if (small_fragments != nb_fragment_num - 1) {
AKANTU_DEBUG_ERROR("The number of small fragments is wrong!");
return EXIT_FAILURE;
}
if (!Math::are_float_equal(total_mass,
small_fragments * frag_theo_mass)) {
AKANTU_DEBUG_ERROR("The mass of small fragments is wrong!");
return EXIT_FAILURE;
}
}
model.dump();
model.dump("cohesive_elements");
/// displace fragments
displaceElements(model, lim, el_size * 2);
model.updateResidual();
lim += el_size;
}
model.dump();
model.dump("cohesive_elements");
/// check centers
- const Array<Real> & fragment_center = model.getFragmentsCenter();
- const Array<Real> & fragment_velocity = model.getFragmentsVelocity();
+ const Array<Real> & fragment_velocity = fragment_manager.getVelocity();
Real initial_position = -L / 2. + el_size / 2.;
for (UInt frag = 0; frag < total_nb_fragment; ++frag) {
Real theoretical_center = initial_position + el_size * frag;
UInt f_index = 0;
while (f_index < total_nb_fragment &&
!Math::are_float_equal(fragment_center(f_index, 0), theoretical_center))
++f_index;
if (f_index == total_nb_fragment) {
AKANTU_DEBUG_ERROR("The fragments' center is wrong!");
return EXIT_FAILURE;
}
f_index = 0;
while (f_index < total_nb_fragment &&
!Math::are_float_equal(fragment_velocity(f_index, 0),
theoretical_center))
++f_index;
if (f_index == total_nb_fragment) {
AKANTU_DEBUG_ERROR("The fragments' velocity is wrong!");
return EXIT_FAILURE;
}
}
finalize();
if (prank == 0)
std::cout << "OK: test_cohesive_buildfragments was passed!" << std::endl;
return EXIT_SUCCESS;
}
void verticalInsertionLimit(SolidMechanicsModelCohesive & model) {
UInt spatial_dimension = model.getSpatialDimension();
const Mesh & mesh_facets = model.getMeshFacets();
const Array<Real> & position = mesh_facets.getNodes();
for (ghost_type_t::iterator gt = ghost_type_t::begin();
gt != ghost_type_t::end();
++gt) {
GhostType ghost_type = *gt;
Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1, ghost_type);
Mesh::type_iterator end = mesh_facets.lastType(spatial_dimension - 1, ghost_type);
for(; it != end; ++it) {
ElementType type = *it;
Array<bool> & check_facets
= model.getElementInserter().getCheckFacets(type, ghost_type);
const Array<UInt> & connectivity = mesh_facets.getConnectivity(type, ghost_type);
UInt nb_nodes_per_facet = connectivity.getNbComponent();
for (UInt f = 0; f < check_facets.getSize(); ++f) {
if (!check_facets(f)) continue;
UInt nb_aligned_nodes = 1;
Real first_node_pos = position(connectivity(f, 0), 0);
for (; nb_aligned_nodes < nb_nodes_per_facet; ++nb_aligned_nodes) {
Real other_node_pos = position(connectivity(f, nb_aligned_nodes), 0);
if (std::abs(first_node_pos - other_node_pos) > 1.e-10)
break;
}
if (nb_aligned_nodes != nb_nodes_per_facet) {
Vector<Real> barycenter(spatial_dimension);
mesh_facets.getBarycenter(f, type, barycenter.storage(), ghost_type);
check_facets(f) = false;
}
}
}
}
}
void displaceElements(SolidMechanicsModelCohesive & model,
const Real lim,
const Real amount) {
UInt spatial_dimension = model.getSpatialDimension();
Array<Real> & displacement = model.getDisplacement();
Mesh & mesh = model.getMesh();
UInt nb_nodes = mesh.getNbNodes();
Array<bool> displaced(nb_nodes);
displaced.clear();
Vector<Real> barycenter(spatial_dimension);
for (ghost_type_t::iterator gt = ghost_type_t::begin();
gt != ghost_type_t::end();
++gt) {
GhostType ghost_type = *gt;
Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type);
Mesh::type_iterator end = mesh.lastType(spatial_dimension, ghost_type);
for(; it != end; ++it) {
ElementType type = *it;
const Array<UInt> & connectivity = mesh.getConnectivity(type, ghost_type);
UInt nb_element = connectivity.getSize();
UInt nb_nodes_per_element = connectivity.getNbComponent();
Array<UInt>::const_iterator<Vector<UInt> > conn_el
= connectivity.begin(nb_nodes_per_element);
for (UInt el = 0; el < nb_element; ++el) {
mesh.getBarycenter(el, type, barycenter.storage(), ghost_type);
if (barycenter(0) < lim) {
const Vector<UInt> & conn = conn_el[el];
for (UInt n = 0; n < nb_nodes_per_element; ++n) {
UInt node = conn(n);
if (!displaced(node)) {
displacement(node, 0) -= amount;
displaced(node) = true;
}
}
}
}
}
}
}