diff --git a/src/fe_engine/element_class.hh b/src/fe_engine/element_class.hh index fce855a3f..5eda6c5ca 100644 --- a/src/fe_engine/element_class.hh +++ b/src/fe_engine/element_class.hh @@ -1,399 +1,399 @@ /** * @file element_class.hh * * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch> * @author Nicolas Richart <nicolas.richart@epfl.ch> * * @date creation: Fri Jun 18 2010 * @date last modification: Thu Jan 21 2016 * * @brief Declaration of the ElementClass main class and the * Integration and Interpolation elements * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see <http://www.gnu.org/licenses/>. * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_types.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_ELEMENT_CLASS_HH__ #define __AKANTU_ELEMENT_CLASS_HH__ namespace akantu { /* -------------------------------------------------------------------------- */ /// default element class structure template <ElementType element_type> struct ElementClassProperty { static const GeometricalType geometrical_type = _gt_not_defined; static const InterpolationType interpolation_type = _itp_not_defined; static const ElementKind element_kind = _ek_regular; static const UInt spatial_dimension = 0; static const GaussIntegrationType gauss_integration_type = _git_not_defined; static const UInt polynomial_degree = 0; }; /// Macro to generate the element class structures for different element types #define AKANTU_DEFINE_ELEMENT_CLASS_PROPERTY(elem_type, geom_type, \ interp_type, elem_kind, sp, \ gauss_int_type, min_int_order) \ template <> struct ElementClassProperty<elem_type> { \ static const GeometricalType geometrical_type = geom_type; \ static const InterpolationType interpolation_type = interp_type; \ static const ElementKind element_kind = elem_kind; \ static const UInt spatial_dimension = sp; \ static const GaussIntegrationType gauss_integration_type = \ gauss_int_type; \ static const UInt polynomial_degree = min_int_order; \ } /* -------------------------------------------------------------------------- */ /* Geometry */ /* -------------------------------------------------------------------------- */ /// Default GeometricalShape structure template <GeometricalType geometrical_type> struct GeometricalShape { static const GeometricalShapeType shape = _gst_point; }; /// Templated GeometricalShape with function contains template <GeometricalShapeType shape> struct GeometricalShapeContains { /// Check if the point (vector in 2 and 3D) at natural coordinate coor template <class vector_type> static inline bool contains(const vector_type & coord); }; /// Macro to generate the GeometricalShape structures for different geometrical /// types #define AKANTU_DEFINE_SHAPE(geom_type, geom_shape) \ template <> struct GeometricalShape<geom_type> { \ static const GeometricalShapeType shape = geom_shape; \ } /* -------------------------------------------------------------------------- */ /// Templated GeometricalElement with function getInradius template <GeometricalType geometrical_type, GeometricalShapeType shape = GeometricalShape<geometrical_type>::shape> class GeometricalElement { public: /// compute the in-radius static inline Real getInradius(__attribute__((unused)) const Matrix<Real> & coord) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// true if the natural coordinates are in the element template <class vector_type> static inline bool contains(const vector_type & coord); public: static AKANTU_GET_MACRO_NOT_CONST(SpatialDimension, spatial_dimension, UInt); static AKANTU_GET_MACRO_NOT_CONST(NbNodesPerElement, nb_nodes_per_element, UInt); static AKANTU_GET_MACRO_NOT_CONST(NbFacetTypes, nb_facet_types, UInt); static inline UInt getNbFacetsPerElement(UInt t); static inline UInt getNbFacetsPerElement(); static inline const MatrixProxy<UInt> getFacetLocalConnectivityPerElement(UInt t = 0); protected: /// Number of nodes per element static UInt nb_nodes_per_element; /// spatial dimension of the element static UInt spatial_dimension; /// number of different facet types static UInt nb_facet_types; /// number of facets for element static UInt nb_facets[]; /// storage of the facet local connectivity static UInt facet_connectivity_vect[]; /// local connectivity of facets static UInt * facet_connectivity[]; private: /// Type of the facet elements static UInt nb_nodes_per_facet[]; }; /* -------------------------------------------------------------------------- */ /* Interpolation */ /* -------------------------------------------------------------------------- */ /// default InterpolationPorperty structure template <InterpolationType interpolation_type> struct InterpolationPorperty { static const InterpolationKind kind = _itk_not_defined; static const UInt nb_nodes_per_element = 0; static const UInt natural_space_dimension = 0; }; /// Macro to generate the InterpolationPorperty structures for different /// interpolation types #define AKANTU_DEFINE_INTERPOLATION_TYPE_PROPERTY(itp_type, itp_kind, \ nb_nodes, ndim) \ template <> struct InterpolationPorperty<itp_type> { \ static const InterpolationKind kind = itp_kind; \ static const UInt nb_nodes_per_element = nb_nodes; \ static const UInt natural_space_dimension = ndim; \ } #include "interpolation_element_tmpl.hh" /* -------------------------------------------------------------------------- */ /// Generic (templated by the enum InterpolationType which specifies the order /// and the dimension of the interpolation) class handling the elemental /// interpolation template <InterpolationType interpolation_type, InterpolationKind kind = InterpolationPorperty<interpolation_type>::kind> class InterpolationElement { public: typedef InterpolationPorperty<interpolation_type> interpolation_property; /// compute the shape values for a given set of points in natural coordinates static inline void computeShapes(const Matrix<Real> & natural_coord, Matrix<Real> & N); /// compute the shape values for a given point in natural coordinates template <class vector_type> static inline void computeShapes(__attribute__((unused)) const vector_type & natural_coord, __attribute__((unused)) vector_type & N) { AKANTU_DEBUG_TO_IMPLEMENT(); } /** * compute @f$ B_{ij} = \frac{\partial N_j}{\partial S_i} @f$ the variation of * shape functions along with variation of natural coordinates on a given set * of points in natural coordinates */ static inline void computeDNDS(const Matrix<Real> & natural_coord, Tensor3<Real> & dnds); /** * compute @f$ B_{ij} = \frac{\partial N_j}{\partial S_i} @f$ the variation of * shape functions along with * variation of natural coordinates on a given point in natural * coordinates */ template <class vector_type, class matrix_type> static inline void computeDNDS(__attribute__((unused)) const vector_type & natural_coord, __attribute__((unused)) matrix_type & dnds) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute jacobian (or integration variable change factor) for a given point /// in the case of spatial_dimension != natural_space_dimension static inline void computeSpecialJacobian(__attribute__((unused)) const Matrix<Real> & J, __attribute__((unused)) Real & jacobians) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// interpolate a field given (arbitrary) natural coordinates static inline void interpolateOnNaturalCoordinates(const Vector<Real> & natural_coords, const Matrix<Real> & nodal_values, Vector<Real> & interpolated); /// interpolate a field given the shape functions on the interpolation point static inline void interpolate(const Matrix<Real> & nodal_values, const Vector<Real> & shapes, Vector<Real> & interpolated); /// interpolate a field given the shape functions on the interpolations points static inline void interpolate(const Matrix<Real> & nodal_values, const Matrix<Real> & shapes, Matrix<Real> & interpolated); /// compute the gradient of a given field on the given natural coordinates static inline void gradientOnNaturalCoordinates(const Vector<Real> & natural_coords, const Matrix<Real> & f, Matrix<Real> & gradient); public: static AKANTU_GET_MACRO_NOT_CONST( ShapeSize, InterpolationPorperty<interpolation_type>::nb_nodes_per_element, UInt); static AKANTU_GET_MACRO_NOT_CONST( ShapeDerivativesSize, (InterpolationPorperty<interpolation_type>::nb_nodes_per_element * InterpolationPorperty<interpolation_type>::natural_space_dimension), UInt); static AKANTU_GET_MACRO_NOT_CONST( NaturalSpaceDimension, InterpolationPorperty<interpolation_type>::natural_space_dimension, UInt); static AKANTU_GET_MACRO_NOT_CONST( NbNodesPerInterpolationElement, InterpolationPorperty<interpolation_type>::nb_nodes_per_element, UInt); }; /* -------------------------------------------------------------------------- */ /* Integration */ /* -------------------------------------------------------------------------- */ template <GaussIntegrationType git_class, UInt nb_points> struct GaussIntegrationTypeData { /// quadrature points in natural coordinates static Real quad_positions[]; /// weights for the Gauss integration static Real quad_weights[]; }; template <ElementType type, UInt n = ElementClassProperty<type>::polynomial_degree> class GaussIntegrationElement { public: static UInt getNbQuadraturePoints(); static const Matrix<Real> getQuadraturePoints(); static const Vector<Real> getWeights(); }; /* -------------------------------------------------------------------------- */ /* ElementClass */ /* -------------------------------------------------------------------------- */ template <ElementType element_type, ElementKind element_kind = ElementClassProperty<element_type>::element_kind> class ElementClass : public GeometricalElement< ElementClassProperty<element_type>::geometrical_type>, public InterpolationElement< ElementClassProperty<element_type>::interpolation_type> { protected: typedef GeometricalElement< ElementClassProperty<element_type>::geometrical_type> geometrical_element; typedef InterpolationElement<ElementClassProperty< element_type>::interpolation_type> interpolation_element; typedef ElementClassProperty<element_type> element_property; typedef typename interpolation_element::interpolation_property interpolation_property; public: /** * compute @f$ J = \frac{\partial x_j}{\partial s_i} @f$ the variation of real * coordinates along with variation of natural coordinates on a given point in * natural coordinates */ static inline void computeJMat(const Matrix<Real> & dnds, const Matrix<Real> & node_coords, Matrix<Real> & J); /** * compute the Jacobian matrix by computing the variation of real coordinates * along with variation of natural coordinates on a given set of points in * natural coordinates */ static inline void computeJMat(const Tensor3<Real> & dnds, const Matrix<Real> & node_coords, Tensor3<Real> & J); /// compute the jacobians of a serie of natural coordinates static inline void computeJacobian(const Matrix<Real> & natural_coords, const Matrix<Real> & node_coords, Vector<Real> & jacobians); /// compute jacobian (or integration variable change factor) for a set of /// points static inline void computeJacobian(const Tensor3<Real> & J, Vector<Real> & jacobians); /// compute jacobian (or integration variable change factor) for a given point static inline void computeJacobian(const Matrix<Real> & J, Real & jacobians); /// compute shape derivatives (input is dxds) for a set of points static inline void computeShapeDerivatives(const Tensor3<Real> & J, const Tensor3<Real> & dnds, Tensor3<Real> & shape_deriv); /// compute shape derivatives (input is dxds) for a given point static inline void computeShapeDerivatives(const Matrix<Real> & J, const Matrix<Real> & dnds, Matrix<Real> & shape_deriv); /// compute the normal of a surface defined by the function f static inline void computeNormalsOnNaturalCoordinates(const Matrix<Real> & coord, Matrix<Real> & f, Matrix<Real> & normals); /// get natural coordinates from real coordinates static inline void inverseMap(const Vector<Real> & real_coords, const Matrix<Real> & node_coords, Vector<Real> & natural_coords, - Real tolerance = 1e-8); + Real tolerance = 1e-10); /// get natural coordinates from real coordinates static inline void inverseMap(const Matrix<Real> & real_coords, const Matrix<Real> & node_coords, Matrix<Real> & natural_coords, - Real tolerance = 1e-8); + Real tolerance = 1e-10); public: static AKANTU_GET_MACRO_NOT_CONST(Kind, element_kind, ElementKind); static AKANTU_GET_MACRO_NOT_CONST( SpatialDimension, ElementClassProperty<element_type>::spatial_dimension, UInt); static AKANTU_GET_MACRO_NOT_CONST(P1ElementType, p1_type, const ElementType &); static const ElementType & getFacetType(UInt t = 0) { return facet_type[t]; } static ElementType * getFacetTypeInternal() { return facet_type; } protected: /// Type of the facet elements static ElementType facet_type[]; /// type of element P1 associated static ElementType p1_type; }; /* -------------------------------------------------------------------------- */ } // akantu #include "element_class_tmpl.hh" /* -------------------------------------------------------------------------- */ namespace akantu { #include "element_class_point_1_inline_impl.cc" #include "element_class_segment_2_inline_impl.cc" #include "element_class_segment_3_inline_impl.cc" #include "element_class_triangle_3_inline_impl.cc" #include "element_class_triangle_6_inline_impl.cc" #include "element_class_tetrahedron_4_inline_impl.cc" #include "element_class_tetrahedron_10_inline_impl.cc" #include "element_class_quadrangle_4_inline_impl.cc" #include "element_class_quadrangle_8_inline_impl.cc" #include "element_class_hexahedron_8_inline_impl.cc" #include "element_class_hexahedron_20_inline_impl.cc" #include "element_class_pentahedron_6_inline_impl.cc" #include "element_class_pentahedron_15_inline_impl.cc" } // akantu /* -------------------------------------------------------------------------- */ #if defined(AKANTU_COHESIVE_ELEMENT) #include "cohesive_element.hh" #endif #if defined(AKANTU_STRUCTURAL_MECHANICS) #include "element_class_structural.hh" #endif #if defined(AKANTU_IGFEM) #include "element_class_igfem.hh" #endif #endif /* __AKANTU_ELEMENT_CLASS_HH__ */ diff --git a/src/fe_engine/element_class_tmpl.hh b/src/fe_engine/element_class_tmpl.hh index e6bdc6db1..51e25a110 100644 --- a/src/fe_engine/element_class_tmpl.hh +++ b/src/fe_engine/element_class_tmpl.hh @@ -1,510 +1,513 @@ /** * @file element_class_tmpl.hh * * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch> * @author Thomas Menouillard <tmenouillard@stucky.ch> * @author Nicolas Richart <nicolas.richart@epfl.ch> * * @date creation: Thu Feb 21 2013 * @date last modification: Thu Oct 22 2015 * * @brief Implementation of the inline templated function of the element class * descriptions * * @section LICENSE * * Copyright (©) 2014, 2015 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 "gauss_integration_tmpl.hh" +#include "element_class.hh" /* -------------------------------------------------------------------------- */ #include <type_traits> /* -------------------------------------------------------------------------- */ +#ifndef __AKANTU_ELEMENT_CLASS_TMPL_HH__ +#define __AKANTU_ELEMENT_CLASS_TMPL_HH__ namespace akantu { - /* -------------------------------------------------------------------------- */ /* GeometricalElement */ /* -------------------------------------------------------------------------- */ template <GeometricalType geometrical_type, GeometricalShapeType shape> inline const MatrixProxy<UInt> GeometricalElement<geometrical_type, shape>::getFacetLocalConnectivityPerElement(UInt t) { return MatrixProxy<UInt>(facet_connectivity[t], nb_facets[t], nb_nodes_per_facet[t]); } /* -------------------------------------------------------------------------- */ template <GeometricalType geometrical_type, GeometricalShapeType shape> inline UInt GeometricalElement<geometrical_type, shape>::getNbFacetsPerElement() { UInt total_nb_facets = 0; for (UInt n = 0; n < nb_facet_types; ++n) { total_nb_facets += nb_facets[n]; } return total_nb_facets; } /* -------------------------------------------------------------------------- */ template <GeometricalType geometrical_type, GeometricalShapeType shape> inline UInt GeometricalElement<geometrical_type, shape>::getNbFacetsPerElement(UInt t) { return nb_facets[t]; } /* -------------------------------------------------------------------------- */ template <GeometricalType geometrical_type, GeometricalShapeType shape> template <class vector_type> inline bool GeometricalElement<geometrical_type, shape>::contains( const vector_type & coords) { return GeometricalShapeContains<shape>::contains(coords); } /* -------------------------------------------------------------------------- */ template <> template <class vector_type> inline bool GeometricalShapeContains<_gst_point>::contains(const vector_type & coords) { return (coords(0) < std::numeric_limits<Real>::epsilon()); } /* -------------------------------------------------------------------------- */ template <> template <class vector_type> inline bool GeometricalShapeContains<_gst_square>::contains(const vector_type & coords) { bool in = true; for (UInt i = 0; i < coords.size() && in; ++i) in &= ((coords(i) >= -(1. + std::numeric_limits<Real>::epsilon())) && (coords(i) <= (1. + std::numeric_limits<Real>::epsilon()))); return in; } /* -------------------------------------------------------------------------- */ template <> template <class vector_type> inline bool GeometricalShapeContains<_gst_triangle>::contains(const vector_type & coords) { bool in = true; Real sum = 0; for (UInt i = 0; (i < coords.size()) && in; ++i) { in &= ((coords(i) >= -(Math::getTolerance())) && (coords(i) <= (1. + Math::getTolerance()))); sum += coords(i); } if (in) return (in && (sum <= (1. + Math::getTolerance()))); return in; } /* -------------------------------------------------------------------------- */ template <> template <class vector_type> inline bool GeometricalShapeContains<_gst_prism>::contains(const vector_type & coords) { bool in = ((coords(0) >= -1.) && (coords(0) <= 1.)); // x in segment [-1, 1] // y and z in triangle in &= ((coords(1) >= 0) && (coords(1) <= 1.)); in &= ((coords(2) >= 0) && (coords(2) <= 1.)); Real sum = coords(1) + coords(2); return (in && (sum <= 1)); } /* -------------------------------------------------------------------------- */ /* InterpolationElement */ /* -------------------------------------------------------------------------- */ template <InterpolationType interpolation_type, InterpolationKind kind> inline void InterpolationElement<interpolation_type, kind>::computeShapes( const Matrix<Real> & natural_coord, Matrix<Real> & N) { UInt nb_points = natural_coord.cols(); for (UInt p = 0; p < nb_points; ++p) { Vector<Real> Np(N(p)); Vector<Real> ncoord_p(natural_coord(p)); computeShapes(ncoord_p, Np); } } /* -------------------------------------------------------------------------- */ template <InterpolationType interpolation_type, InterpolationKind kind> inline void InterpolationElement<interpolation_type, kind>::computeDNDS( const Matrix<Real> & natural_coord, Tensor3<Real> & dnds) { UInt nb_points = natural_coord.cols(); for (UInt p = 0; p < nb_points; ++p) { Matrix<Real> dnds_p(dnds(p)); Vector<Real> ncoord_p(natural_coord(p)); computeDNDS(ncoord_p, dnds_p); } } /* -------------------------------------------------------------------------- */ /** * interpolate on a point a field for which values are given on the * node of the element using the shape functions at this interpolation point * * @param nodal_values values of the function per node @f$ f_{ij} = f_{n_i j} *@f$ so it should be a matrix of size nb_nodes_per_element @f$\times@f$ *nb_degree_of_freedom * @param shapes value of shape functions at the interpolation point * @param interpolated interpolated value of f @f$ f_j(\xi) = \sum_i f_{n_i j} *N_i @f$ */ template <InterpolationType interpolation_type, InterpolationKind kind> inline void InterpolationElement<interpolation_type, kind>::interpolate( const Matrix<Real> & nodal_values, const Vector<Real> & shapes, Vector<Real> & interpolated) { Matrix<Real> interpm(interpolated.storage(), nodal_values.rows(), 1); Matrix<Real> shapesm( shapes.storage(), InterpolationPorperty<interpolation_type>::nb_nodes_per_element, 1); interpm.mul<false, false>(nodal_values, shapesm); } /* -------------------------------------------------------------------------- */ /** * interpolate on several points a field for which values are given on the * node of the element using the shape functions at the interpolation point * * @param nodal_values values of the function per node @f$ f_{ij} = f_{n_i j} *@f$ so it should be a matrix of size nb_nodes_per_element @f$\times@f$ *nb_degree_of_freedom * @param shapes value of shape functions at the interpolation point * @param interpolated interpolated values of f @f$ f_j(\xi) = \sum_i f_{n_i j} *N_i @f$ */ template <InterpolationType interpolation_type, InterpolationKind kind> inline void InterpolationElement<interpolation_type, kind>::interpolate( const Matrix<Real> & nodal_values, const Matrix<Real> & shapes, Matrix<Real> & interpolated) { UInt nb_points = shapes.cols(); for (UInt p = 0; p < nb_points; ++p) { Vector<Real> Np(shapes(p)); Vector<Real> interpolated_p(interpolated(p)); interpolate(nodal_values, Np, interpolated_p); } } /* -------------------------------------------------------------------------- */ /** * interpolate the field on a point given in natural coordinates the field which * values are given on the node of the element * * @param natural_coords natural coordinates of point where to interpolate \xi * @param nodal_values values of the function per node @f$ f_{ij} = f_{n_i j} *@f$ so it should be a matrix of size nb_nodes_per_element @f$\times@f$ *nb_degree_of_freedom * @param interpolated interpolated value of f @f$ f_j(\xi) = \sum_i f_{n_i j} *N_i @f$ */ template <InterpolationType interpolation_type, InterpolationKind kind> inline void InterpolationElement<interpolation_type, kind>::interpolateOnNaturalCoordinates( const Vector<Real> & natural_coords, const Matrix<Real> & nodal_values, Vector<Real> & interpolated) { Vector<Real> shapes( InterpolationPorperty<interpolation_type>::nb_nodes_per_element); computeShapes(natural_coords, shapes); interpolate(nodal_values, shapes, interpolated); } /* -------------------------------------------------------------------------- */ /// @f$ gradient_{ij} = \frac{\partial f_j}{\partial s_i} = \sum_k /// \frac{\partial N_k}{\partial s_i}f_{j n_k} @f$ template <InterpolationType interpolation_type, InterpolationKind kind> inline void InterpolationElement<interpolation_type, kind>::gradientOnNaturalCoordinates( const Vector<Real> & natural_coords, const Matrix<Real> & f, Matrix<Real> & gradient) { Matrix<Real> dnds( InterpolationPorperty<interpolation_type>::natural_space_dimension, InterpolationPorperty<interpolation_type>::nb_nodes_per_element); computeDNDS(natural_coords, dnds); gradient.mul<false, true>(f, dnds); } /* -------------------------------------------------------------------------- */ /* ElementClass */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template <ElementType type, ElementKind kind> inline void ElementClass<type, kind>::computeJMat(const Tensor3<Real> & dnds, const Matrix<Real> & node_coords, Tensor3<Real> & J) { UInt nb_points = dnds.size(2); for (UInt p = 0; p < nb_points; ++p) { Matrix<Real> J_p(J(p)); Matrix<Real> dnds_p(dnds(p)); computeJMat(dnds_p, node_coords, J_p); } } /* -------------------------------------------------------------------------- */ template <ElementType type, ElementKind kind> inline void ElementClass<type, kind>::computeJMat(const Matrix<Real> & dnds, const Matrix<Real> & node_coords, Matrix<Real> & J) { /// @f$ J = dxds = dnds * x @f$ J.mul<false, true>(dnds, node_coords); } /* -------------------------------------------------------------------------- */ template <ElementType type, ElementKind kind> inline void ElementClass<type, kind>::computeJacobian(const Matrix<Real> & natural_coords, const Matrix<Real> & node_coords, Vector<Real> & jacobians) { UInt nb_points = natural_coords.cols(); Matrix<Real> dnds(interpolation_property::natural_space_dimension, interpolation_property::nb_nodes_per_element); Matrix<Real> J(natural_coords.rows(), node_coords.rows()); for (UInt p = 0; p < nb_points; ++p) { Vector<Real> ncoord_p(natural_coords(p)); interpolation_element::computeDNDS(ncoord_p, dnds); computeJMat(dnds, node_coords, J); computeJacobian(J, jacobians(p)); } } /* -------------------------------------------------------------------------- */ template <ElementType type, ElementKind kind> inline void ElementClass<type, kind>::computeJacobian(const Tensor3<Real> & J, Vector<Real> & jacobians) { UInt nb_points = J.size(2); for (UInt p = 0; p < nb_points; ++p) { computeJacobian(J(p), jacobians(p)); } } /* -------------------------------------------------------------------------- */ template <ElementType type, ElementKind kind> inline void ElementClass<type, kind>::computeJacobian(const Matrix<Real> & J, Real & jacobians) { if (J.rows() == J.cols()) { jacobians = Math::det<element_property::spatial_dimension>(J.storage()); } else { interpolation_element::computeSpecialJacobian(J, jacobians); } } /* -------------------------------------------------------------------------- */ template <ElementType type, ElementKind kind> inline void ElementClass<type, kind>::computeShapeDerivatives(const Tensor3<Real> & J, const Tensor3<Real> & dnds, Tensor3<Real> & shape_deriv) { UInt nb_points = J.size(2); for (UInt p = 0; p < nb_points; ++p) { Matrix<Real> shape_deriv_p(shape_deriv(p)); computeShapeDerivatives(J(p), dnds(p), shape_deriv_p); } } /* -------------------------------------------------------------------------- */ template <ElementType type, ElementKind kind> inline void ElementClass<type, kind>::computeShapeDerivatives(const Matrix<Real> & J, const Matrix<Real> & dnds, Matrix<Real> & shape_deriv) { Matrix<Real> inv_J(J.rows(), J.cols()); Math::inv<element_property::spatial_dimension>(J.storage(), inv_J.storage()); shape_deriv.mul<false, false>(inv_J, dnds); } /* -------------------------------------------------------------------------- */ template <ElementType type, ElementKind kind> inline void ElementClass<type, kind>::computeNormalsOnNaturalCoordinates( const Matrix<Real> & coord, Matrix<Real> & f, Matrix<Real> & normals) { UInt dimension = normals.rows(); UInt nb_points = coord.cols(); AKANTU_DEBUG_ASSERT((dimension - 1) == interpolation_property::natural_space_dimension, "cannot extract a normal because of dimension mismatch " << dimension - 1 << " " << interpolation_property::natural_space_dimension); Matrix<Real> J(dimension, interpolation_property::natural_space_dimension); for (UInt p = 0; p < nb_points; ++p) { interpolation_element::gradientOnNaturalCoordinates(coord(p), f, J); if (dimension == 2) { Math::normal2(J.storage(), normals(p).storage()); } if (dimension == 3) { Math::normal3(J(0).storage(), J(1).storage(), normals(p).storage()); } } } /* ------------------------------------------------------------------------- */ /** * In the non linear cases we need to iterate to find the natural coordinates *@f$\xi@f$ * provided real coordinates @f$x@f$. * * We want to solve: @f$ x- \phi(\xi) = 0@f$ with @f$\phi(\xi) = \sum_I N_I(\xi) *x_I@f$ * the mapping function which uses the nodal coordinates @f$x_I@f$. * * To that end we use the Newton method and the following series: * * @f$ \frac{\partial \phi(x_k)}{\partial \xi} \left( \xi_{k+1} - \xi_k \right) *= x - \phi(x_k)@f$ * * When we consider elements embedded in a dimension higher than them (2D *triangle in a 3D space for example) * @f$ J = \frac{\partial \phi(\xi_k)}{\partial \xi}@f$ is of dimension *@f$dim_{space} \times dim_{elem}@f$ which * is not invertible in most cases. Rather we can solve the problem: * * @f$ J^T J \left( \xi_{k+1} - \xi_k \right) = J^T \left( x - \phi(\xi_k) *\right) @f$ * * So that * * @f$ d\xi = \xi_{k+1} - \xi_k = (J^T J)^{-1} J^T \left( x - \phi(\xi_k) *\right) @f$ * * So that if the series converges we have: * * @f$ 0 = J^T \left( \phi(\xi_\infty) - x \right) @f$ * * And we see that this is ill-posed only if @f$ J^T x = 0@f$ which means that *the vector provided * is normal to any tangent which means it is outside of the element itself. * * @param real_coords: the real coordinates the natural coordinates are sought *for * @param node_coords: the coordinates of the nodes forming the element * @param natural_coords: output->the sought natural coordinates * @param spatial_dimension: spatial dimension of the problem * **/ template <ElementType type, ElementKind kind> inline void ElementClass<type, kind>::inverseMap( const Vector<Real> & real_coords, const Matrix<Real> & node_coords, Vector<Real> & natural_coords, Real tolerance) { UInt spatial_dimension = real_coords.size(); UInt dimension = natural_coords.size(); // matrix copy of the real_coords Matrix<Real> mreal_coords(real_coords.storage(), spatial_dimension, 1); // initial guess // Matrix<Real> natural_guess(natural_coords.storage(), dimension, 1); natural_coords.clear(); // real space coordinates provided by initial guess Matrix<Real> physical_guess(dimension, 1); // objective function f = real_coords - physical_guess Matrix<Real> f(dimension, 1); // dnds computed on the natural_guess // Matrix<Real> dnds(interpolation_element::nb_nodes_per_element, // spatial_dimension); // J Jacobian matrix computed on the natural_guess Matrix<Real> J(spatial_dimension, dimension); // G = J^t * J Matrix<Real> G(spatial_dimension, spatial_dimension); // Ginv = G^{-1} Matrix<Real> Ginv(spatial_dimension, spatial_dimension); // J = Ginv * J^t Matrix<Real> F(spatial_dimension, dimension); // dxi = \xi_{k+1} - \xi in the iterative process Matrix<Real> dxi(spatial_dimension, 1); /* --------------------------- */ /* init before iteration loop */ /* --------------------------- */ // do interpolation - Vector<Real> physical_guess_v(physical_guess.storage(), dimension); - interpolation_element::interpolateOnNaturalCoordinates( - natural_coords, node_coords, physical_guess_v); + auto update_f = [&f, &physical_guess, &natural_coords, &node_coords, &mreal_coords, + dimension]() { + Vector<Real> physical_guess_v(physical_guess.storage(), dimension); + interpolation_element::interpolateOnNaturalCoordinates( + natural_coords, node_coords, physical_guess_v); + + // compute initial objective function value f = real_coords - physical_guess + f = mreal_coords; + f -= physical_guess; - // compute initial objective function value f = real_coords - physical_guess - f = mreal_coords; - f -= physical_guess; + // compute initial error + auto error = f.norm<L_2>(); + return error; + }; - // compute initial error - Real inverse_map_error = f.norm<L_2>(); + auto inverse_map_error = update_f(); /* --------------------------- */ /* iteration loop */ /* --------------------------- */ while (tolerance < inverse_map_error) { // compute J^t interpolation_element::gradientOnNaturalCoordinates(natural_coords, node_coords, J); // compute G G.mul<true, false>(J, J); // inverse G Ginv.inverse(G); // compute F F.mul<false, true>(Ginv, J); // compute increment dxi.mul<false, false>(F, f); // update our guess natural_coords += Vector<Real>(dxi(0)); - // interpolate - interpolation_element::interpolateOnNaturalCoordinates( - natural_coords, node_coords, physical_guess_v); - - // compute error - f = mreal_coords; - f -= physical_guess; - inverse_map_error = f.norm<L_2>(); + inverse_map_error = update_f(); } // memcpy(natural_coords.storage(), natural_guess.storage(), sizeof(Real) * // natural_coords.size()); } /* -------------------------------------------------------------------------- */ template <ElementType type, ElementKind kind> inline void ElementClass<type, kind>::inverseMap( const Matrix<Real> & real_coords, const Matrix<Real> & node_coords, Matrix<Real> & natural_coords, Real tolerance) { UInt nb_points = real_coords.cols(); for (UInt p = 0; p < nb_points; ++p) { Vector<Real> X(real_coords(p)); Vector<Real> ncoord_p(natural_coords(p)); inverseMap(X, node_coords, ncoord_p, tolerance); } } -} // akantu +} // namespace akantu + +#endif /* __AKANTU_ELEMENT_CLASS_TMPL_HH__ */ diff --git a/src/fe_engine/fe_engine_template.hh b/src/fe_engine/fe_engine_template.hh index dc8a0c383..fcca5d08e 100644 --- a/src/fe_engine/fe_engine_template.hh +++ b/src/fe_engine/fe_engine_template.hh @@ -1,385 +1,394 @@ /** * @file fe_engine_template.hh * * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> * @author Nicolas Richart <nicolas.richart@epfl.ch> * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Dec 08 2015 * * @brief templated class that calls integration and shape objects * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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_FE_ENGINE_TEMPLATE_HH__ #define __AKANTU_FE_ENGINE_TEMPLATE_HH__ /* -------------------------------------------------------------------------- */ #include "fe_engine.hh" #include "integrator.hh" #include "shape_functions.hh" /* -------------------------------------------------------------------------- */ #include <type_traits> /* -------------------------------------------------------------------------- */ namespace akantu { class DOFManager; namespace fe_engine { namespace details { template <ElementKind> struct AssembleLumpedTemplateHelper; template <ElementKind> struct AssembleFieldMatrixHelper; } // namespace details } // namespace fe_engine template <ElementKind, typename> struct AssembleFieldMatrixStructHelper; struct DefaultIntegrationOrderFunctor { template <ElementType type> static inline constexpr int getOrder() { return ElementClassProperty<type>::polynomial_degree; } }; /* -------------------------------------------------------------------------- */ template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind = _ek_regular, class IntegrationOrderFunctor = DefaultIntegrationOrderFunctor> class FEEngineTemplate : public FEEngine { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef I<kind, IntegrationOrderFunctor> Integ; typedef S<kind> Shape; FEEngineTemplate(Mesh & mesh, UInt spatial_dimension = _all_dimensions, ID id = "fem", MemoryID memory_id = 0); virtual ~FEEngineTemplate(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// pre-compute all the shape functions, their derivatives and the jacobians void initShapeFunctions(const GhostType & ghost_type = _not_ghost) override; 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 override; /// 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 override; /// 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 override; /// integrate partially around an integration point (@f$ intf_q = f_q * J_q * /// w_q @f$) void integrateOnIntegrationPoints( 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 override; /// interpolate on a phyiscal point inside an element void interpolate(const Vector<Real> & real_coords, const Matrix<Real> & nodal_values, Vector<Real> & interpolated, const Element & element) const override; /// get the number of integration points UInt getNbIntegrationPoints( const ElementType & type, const GhostType & ghost_type = _not_ghost) const override; /// get shapes precomputed const Array<Real> & getShapes(const ElementType & type, const GhostType & ghost_type = _not_ghost, UInt id = 0) const override; /// get the derivatives of shapes const Array<Real> & getShapesDerivatives(const ElementType & type, const GhostType & ghost_type = _not_ghost, UInt id = 0) const override; /// get integration points const inline Matrix<Real> & getIntegrationPoints( const ElementType & type, const GhostType & ghost_type = _not_ghost) const override; /* ------------------------------------------------------------------------ */ /* Shape method bridges */ /* ------------------------------------------------------------------------ */ /// compute the gradient of a nodal field on the integration points void gradientOnIntegrationPoints( 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 override; /// interpolate a nodal field on the integration points void interpolateOnIntegrationPoints( 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 override; /// interpolate a nodal field on the integration points given a /// by_element_type void interpolateOnIntegrationPoints( const Array<Real> & u, ElementTypeMapArray<Real> & uq, const ElementTypeMapArray<UInt> * filter_elements = NULL) const override; /// compute the position of integration points given by an element_type_map /// from nodes position inline void computeIntegrationPointsCoordinates( ElementTypeMapArray<Real> & quadrature_points_coordinates, const ElementTypeMapArray<UInt> * filter_elements = NULL) const override; /// compute the position of integration points from nodes position inline void computeIntegrationPointsCoordinates( Array<Real> & quadrature_points_coordinates, const ElementType & type, const GhostType & ghost_type = _not_ghost, const Array<UInt> & filter_elements = empty_filter) const override; /// interpolate field at given position (interpolation_points) from given /// values of this field at integration points (field) inline void interpolateElementalFieldFromIntegrationPoints( const ElementTypeMapArray<Real> & field, const ElementTypeMapArray<Real> & interpolation_points_coordinates, ElementTypeMapArray<Real> & result, const GhostType ghost_type, const ElementTypeMapArray<UInt> * element_filter) const override; /// Interpolate field at given position from given values of this field at /// integration points (field) /// using matrices precomputed with /// initElementalFieldInterplationFromIntegrationPoints inline void interpolateElementalFieldFromIntegrationPoints( const ElementTypeMapArray<Real> & field, const ElementTypeMapArray<Real> & interpolation_points_coordinates_matrices, const ElementTypeMapArray<Real> & quad_points_coordinates_inv_matrices, ElementTypeMapArray<Real> & result, const GhostType ghost_type, const ElementTypeMapArray<UInt> * element_filter) const override; /// Build pre-computed matrices for interpolation of field form integration /// points at other given positions (interpolation_points) inline void initElementalFieldInterpolationFromIntegrationPoints( const ElementTypeMapArray<Real> & interpolation_points_coordinates, ElementTypeMapArray<Real> & interpolation_points_coordinates_matrices, ElementTypeMapArray<Real> & quad_points_coordinates_inv_matrices, const ElementTypeMapArray<UInt> * element_filter = NULL) const override; /// 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 override; /// compute the shape derivatives on a provided point inline void computeShapeDerivatives( const Vector<Real> & real__coords, UInt element, const ElementType & type, Matrix<Real> & shape_derivatives, const GhostType & ghost_type = _not_ghost) const override; /* ------------------------------------------------------------------------ */ /* Other methods */ /* ------------------------------------------------------------------------ */ /// pre-compute normals on integration points void computeNormalsOnIntegrationPoints( const GhostType & ghost_type = _not_ghost) override; void computeNormalsOnIntegrationPoints( const Array<Real> & field, const GhostType & ghost_type = _not_ghost) override; void computeNormalsOnIntegrationPoints( const Array<Real> & field, Array<Real> & normal, const ElementType & type, const GhostType & ghost_type = _not_ghost) const override; template <ElementType type> void computeNormalsOnIntegrationPoints(const Array<Real> & field, Array<Real> & normal, const GhostType & ghost_type) const; +private: + // To avoid a weird full specialization of a method in a non specalized class + void + computeNormalsOnIntegrationPointsPoint1(const Array<Real> &, + Array<Real> & normal, + const GhostType & ghost_type) const; + +public: /// function to print the contain of the class void printself(std::ostream & stream, int indent = 0) const override; /// assemble a field as a lumped matrix (ex. rho in lumped mass) template <class Functor> void assembleFieldLumped(const Functor & field_funct, const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager, ElementType type, const GhostType & ghost_type) const; /// assemble a field as a matrix (ex. rho to mass matrix) template <class Functor> void assembleFieldMatrix(const Functor & field_funct, const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager, ElementType type, const GhostType & ghost_type) const; #ifdef AKANTU_STRUCTURAL_MECHANICS /// assemble a field as a matrix (ex. rho to mass matrix) void assembleFieldMatrix(const Array<Real> & field_1, UInt nb_degree_of_freedom, SparseMatrix & M, Array<Real> * n, ElementTypeMapArray<Real> & rotation_mat, const ElementType & type, const GhostType & ghost_type = _not_ghost) const; /// compute shapes function in a matrix for structural elements void computeShapesMatrix(const ElementType & type, UInt nb_degree_of_freedom, UInt nb_nodes_per_element, Array<Real> * n, UInt id, UInt degree_to_interpolate, UInt degree_interpolated, const bool sign, const GhostType & ghost_type = _not_ghost) const; #endif private: friend struct fe_engine::details::AssembleLumpedTemplateHelper<kind>; friend struct fe_engine::details::AssembleFieldMatrixHelper<kind>; friend struct AssembleFieldMatrixStructHelper<kind, void>; /// templated function to compute the scaling to assemble a lumped matrix template <class Functor, ElementType type> void assembleFieldLumped(const Functor & field_funct, const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager, 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, const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager, 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, const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager, const GhostType & ghost_type) const; /// assemble a field as a matrix (ex. rho to mass matrix) template <class Functor, ElementType type> void assembleFieldMatrix(const Functor & field_funct, const ID & matrix_id, const ID & dof_id, DOFManager & dof_manager, const GhostType & ghost_type) const; #ifdef AKANTU_STRUCTURAL_MECHANICS /// assemble a field as a matrix for structural elements (ex. rho to mass /// matrix) template <ElementType type> void assembleFieldMatrix(const Array<Real> & field_1, UInt nb_degree_of_freedom, SparseMatrix & M, Array<Real> * n, ElementTypeMapArray<Real> & rotation_mat, __attribute__((unused)) const GhostType & ghost_type) const; #endif /* ------------------------------------------------------------------------ */ /* Mesh Event Handler interface */ /* ------------------------------------------------------------------------ */ public: void onElementsAdded(const Array<Element> &, const NewElementsEvent &) override; void onElementsRemoved(const Array<Element> &, const ElementTypeMapArray<UInt> &, const RemovedElementsEvent &) override; void onElementsChanged(const Array<Element> &, const Array<Element> &, const ElementTypeMapArray<UInt> &, const ChangedElementsEvent &) override; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the shape class (probably useless: see getShapeFunction) const ShapeFunctions & getShapeFunctionsInterface() const override { return shape_functions; }; /// get the shape class const Shape & getShapeFunctions() const { return shape_functions; }; /// get the integrator class (probably useless: see getIntegrator) const Integrator & getIntegratorInterface() const override { return integrator; }; /// get the integrator class const Integ & getIntegrator() const { return integrator; }; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: Integ integrator; Shape shape_functions; }; } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "fe_engine_template_tmpl.hh" #include "fe_engine_template_tmpl_field.hh" /* -------------------------------------------------------------------------- */ /* Shape Linked specialization */ /* -------------------------------------------------------------------------- */ #if defined(AKANTU_STRUCTURAL_MECHANICS) #include "fe_engine_template_tmpl_struct.hh" #endif /* -------------------------------------------------------------------------- */ /* Shape IGFEM specialization */ /* -------------------------------------------------------------------------- */ #if defined(AKANTU_IGFEM) #include "fe_engine_template_tmpl_igfem.hh" #endif #endif /* __AKANTU_FE_ENGINE_TEMPLATE_HH__ */ diff --git a/src/fe_engine/fe_engine_template_cohesive.cc b/src/fe_engine/fe_engine_template_cohesive.cc index 0e5616dee..8e98c9058 100644 --- a/src/fe_engine/fe_engine_template_cohesive.cc +++ b/src/fe_engine/fe_engine_template_cohesive.cc @@ -1,175 +1,176 @@ /** * @file fe_engine_template_cohesive.cc * * @author Nicolas Richart <nicolas.richart@epfl.ch> * @author Marco Vocialta <marco.vocialta@epfl.ch> * * @date creation: Wed Oct 31 2012 * @date last modification: Thu Oct 15 2015 * * @brief Specialization of the FEEngineTemplate for cohesive element * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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 "fe_engine_template.hh" #include "shape_cohesive.hh" +#include "integrator_gauss.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /* compatibility functions */ /* -------------------------------------------------------------------------- */ template <> Real FEEngineTemplate<IntegratorGauss, ShapeLagrange, _ek_cohesive, DefaultIntegrationOrderFunctor>:: integrate(const Array<Real> & f, const ElementType & type, const GhostType & ghost_type, const Array<UInt> & filter_elements) const { AKANTU_DEBUG_IN(); #ifndef AKANTU_NDEBUG UInt nb_element = mesh.getNbElement(type, ghost_type); if (filter_elements != empty_filter) nb_element = filter_elements.getSize(); UInt nb_quadrature_points = getNbIntegrationPoints(type); AKANTU_DEBUG_ASSERT(f.getSize() == nb_element * nb_quadrature_points, "The vector f(" << f.getID() << ") has not the good size."); 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.integrate<type>(f, ghost_type, filter_elements); AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(INTEGRATE); #undef INTEGRATE AKANTU_DEBUG_OUT(); return integral; } /* -------------------------------------------------------------------------- */ template <> void FEEngineTemplate<IntegratorGauss, ShapeLagrange, _ek_cohesive, DefaultIntegrationOrderFunctor>:: 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 { #ifndef AKANTU_NDEBUG UInt nb_element = mesh.getNbElement(type, ghost_type); if (filter_elements == filter_elements) nb_element = filter_elements.getSize(); UInt nb_quadrature_points = getNbIntegrationPoints(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."); AKANTU_DEBUG_ASSERT(intf.getSize() == nb_element, "The vector intf(" << intf.getID() << ") has not the good size."); #endif #define INTEGRATE(type) \ integrator.integrate<type>(f, intf, nb_degree_of_freedom, ghost_type, \ filter_elements); AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(INTEGRATE); #undef INTEGRATE } /* -------------------------------------------------------------------------- */ template <> void FEEngineTemplate<IntegratorGauss, ShapeLagrange, _ek_cohesive, DefaultIntegrationOrderFunctor>:: gradientOnIntegrationPoints(const Array<Real> &, Array<Real> &, const UInt, const ElementType &, const GhostType &, const Array<UInt> &) const { AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ template <> template <> void FEEngineTemplate<IntegratorGauss, ShapeLagrange, _ek_cohesive, DefaultIntegrationOrderFunctor>:: computeNormalsOnIntegrationPoints<_cohesive_1d_2>( const Array<Real> &, 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 1D cohesive elements!"); const ElementType type = _cohesive_1d_2; const ElementType facet_type = Mesh::getFacetType(type); UInt nb_element = mesh.getConnectivity(type, ghost_type).getSize(); normal.resize(nb_element); Array<Element> & facets = mesh.getMeshFacets().getSubelementToElement(type, ghost_type); Array<std::vector<Element> > & segments = mesh.getMeshFacets().getElementToSubelement(facet_type, ghost_type); Real values[2]; for (UInt elem = 0; elem < nb_element; ++elem) { for (UInt p = 0; p < 2; ++p) { Element f = facets(elem, p); Element seg = segments(f.element)[0]; mesh.getBarycenter(seg.element, seg.type, &(values[p]), seg.ghost_type); } Real difference = values[0] - values[1]; AKANTU_DEBUG_ASSERT(difference != 0., "Error in normal computation for cohesive elements"); normal(elem) = difference / std::abs(difference); } AKANTU_DEBUG_OUT(); } } // akantu diff --git a/src/fe_engine/fe_engine_template_tmpl.hh b/src/fe_engine/fe_engine_template_tmpl.hh index da5a03be7..2f95e57c5 100644 --- a/src/fe_engine/fe_engine_template_tmpl.hh +++ b/src/fe_engine/fe_engine_template_tmpl.hh @@ -1,1306 +1,1301 @@ /** * @file fe_engine_template_tmpl.hh * * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> * @author Mauro Corrado <mauro.corrado@epfl.ch> * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch> * @author Nicolas Richart <nicolas.richart@epfl.ch> * @author Marco Vocialta <marco.vocialta@epfl.ch> * * @date creation: Tue Feb 15 2011 * @date last modification: Thu Nov 19 2015 * * @brief Template implementation of FEEngineTemplate * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see <http://www.gnu.org/licenses/>. * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "dof_manager.hh" #include "fe_engine_template.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::FEEngineTemplate( Mesh & mesh, UInt spatial_dimension, ID id, MemoryID memory_id) : FEEngine(mesh, spatial_dimension, id, memory_id), integrator(mesh, id, memory_id), shape_functions(mesh, id, memory_id) {} /* -------------------------------------------------------------------------- */ template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::~FEEngineTemplate() {} /* -------------------------------------------------------------------------- */ /** * Helper class to be able to write a partial specialization on the element kind */ namespace fe_engine { namespace details { template <ElementKind kind> struct GradientOnIntegrationPointsHelper { template <class S> static void call(const S &, Mesh &, const Array<Real> &, Array<Real> &, const UInt, const ElementType &, const GhostType &, const Array<UInt> &) { AKANTU_DEBUG_TO_IMPLEMENT(); } }; #define COMPUTE_GRADIENT(type) \ if (element_dimension == ElementClass<type>::getSpatialDimension()) \ shape_functions.template gradientOnIntegrationPoints<type>( \ u, nablauq, nb_degree_of_freedom, ghost_type, filter_elements); #define AKANTU_SPECIALIZE_GRADIENT_ON_INTEGRATION_POINTS_HELPER(kind) \ template <> struct GradientOnIntegrationPointsHelper<kind> { \ template <class S> \ static void call(const S & shape_functions, Mesh & mesh, \ 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) { \ UInt element_dimension = mesh.getSpatialDimension(type); \ AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_GRADIENT, kind); \ } \ }; AKANTU_BOOST_ALL_KIND_LIST( AKANTU_SPECIALIZE_GRADIENT_ON_INTEGRATION_POINTS_HELPER, AKANTU_FE_ENGINE_LIST_GRADIENT_ON_INTEGRATION_POINTS) #undef AKANTU_SPECIALIZE_GRADIENT_ON_INTEGRATION_POINTS_HELPER #undef COMPUTE_GRADIENT } // namespace details } // namespace fe_engine template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>:: gradientOnIntegrationPoints(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.getIntegrationPoints(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."); #endif nablauq.resize(nb_element * nb_points); fe_engine::details::GradientOnIntegrationPointsHelper<kind>::call( shape_functions, mesh, u, nablauq, nb_degree_of_freedom, type, ghost_type, filter_elements); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::initShapeFunctions( const GhostType & ghost_type) { initShapeFunctions(mesh.getNodes(), ghost_type); } /* -------------------------------------------------------------------------- */ template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::initShapeFunctions( const Array<Real> & nodes, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); for (auto & type : mesh.elementTypes(element_dimension, ghost_type, kind)) { integrator.initIntegrator(nodes, type, ghost_type); const auto & control_points = getIntegrationPoints(type, ghost_type); shape_functions.initShapeFunctions(nodes, control_points, type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Helper class to be able to write a partial specialization on the element kind */ namespace fe_engine { namespace details { template <ElementKind kind> struct IntegrateHelper {}; #define INTEGRATE(type) \ integrator.template integrate<type>(f, intf, nb_degree_of_freedom, \ ghost_type, filter_elements); #define AKANTU_SPECIALIZE_INTEGRATE_HELPER(kind) \ template <> struct IntegrateHelper<kind> { \ template <class I> \ static void call(const I & integrator, const Array<Real> & f, \ Array<Real> & intf, UInt nb_degree_of_freedom, \ const ElementType & type, const GhostType & ghost_type, \ const Array<UInt> & filter_elements) { \ AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTEGRATE, kind); \ } \ }; AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_INTEGRATE_HELPER) #undef AKANTU_SPECIALIZE_INTEGRATE_HELPER #undef INTEGRATE } // namespace details } // namespace fe_engine template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::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 = getNbIntegrationPoints(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); fe_engine::details::IntegrateHelper<kind>::call(integrator, f, intf, nb_degree_of_freedom, type, ghost_type, filter_elements); } /* -------------------------------------------------------------------------- */ /** * Helper class to be able to write a partial specialization on the element kind */ namespace fe_engine { namespace details { template <ElementKind kind> struct IntegrateScalarHelper {}; #define INTEGRATE(type) \ integral = \ integrator.template integrate<type>(f, ghost_type, filter_elements); #define AKANTU_SPECIALIZE_INTEGRATE_SCALAR_HELPER(kind) \ template <> struct IntegrateScalarHelper<kind> { \ template <class I> \ static Real call(const I & integrator, const Array<Real> & f, \ const ElementType & type, const GhostType & ghost_type, \ const Array<UInt> & filter_elements) { \ Real integral = 0.; \ AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTEGRATE, kind); \ return integral; \ } \ }; AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_INTEGRATE_SCALAR_HELPER) #undef AKANTU_SPECIALIZE_INTEGRATE_SCALAR_HELPER #undef INTEGRATE } // namespace details } // namespace fe_engine template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> Real FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::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 = getNbIntegrationPoints(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 = fe_engine::details::IntegrateScalarHelper<kind>::call( integrator, f, type, ghost_type, filter_elements); AKANTU_DEBUG_OUT(); return integral; } /* -------------------------------------------------------------------------- */ /** * Helper class to be able to write a partial specialization on the element kind */ namespace fe_engine { namespace details { template <ElementKind kind> struct IntegrateScalarOnOneElementHelper {}; #define INTEGRATE(type) \ res = integrator.template integrate<type>(f, index, ghost_type); #define AKANTU_SPECIALIZE_INTEGRATE_SCALAR_ON_ONE_ELEMENT_HELPER(kind) \ template <> struct IntegrateScalarOnOneElementHelper<kind> { \ template <class I> \ static Real call(const I & integrator, const Vector<Real> & f, \ const ElementType & type, UInt index, \ const GhostType & ghost_type) { \ Real res = 0.; \ AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTEGRATE, kind); \ return res; \ } \ }; AKANTU_BOOST_ALL_KIND( AKANTU_SPECIALIZE_INTEGRATE_SCALAR_ON_ONE_ELEMENT_HELPER) #undef AKANTU_SPECIALIZE_INTEGRATE_SCALAR_ON_ONE_ELEMENT_HELPER #undef INTEGRATE } // namespace details } // namespace fe_engine template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> Real FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::integrate( const Vector<Real> & f, const ElementType & type, UInt index, const GhostType & ghost_type) const { Real res = fe_engine::details::IntegrateScalarOnOneElementHelper<kind>::call( integrator, f, type, index, ghost_type); return res; } /* -------------------------------------------------------------------------- */ /** * Helper class to be able to write a partial specialization on the element kind */ namespace fe_engine { namespace details { template <ElementKind kind> struct IntegrateOnIntegrationPointsHelper {}; #define INTEGRATE(type) \ integrator.template integrateOnIntegrationPoints<type>( \ f, intf, nb_degree_of_freedom, ghost_type, filter_elements); #define AKANTU_SPECIALIZE_INTEGRATE_ON_INTEGRATION_POINTS_HELPER(kind) \ template <> struct IntegrateOnIntegrationPointsHelper<kind> { \ template <class I> \ static void call(const I & integrator, const Array<Real> & f, \ Array<Real> & intf, UInt nb_degree_of_freedom, \ const ElementType & type, const GhostType & ghost_type, \ const Array<UInt> & filter_elements) { \ AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTEGRATE, kind); \ } \ }; AKANTU_BOOST_ALL_KIND( AKANTU_SPECIALIZE_INTEGRATE_ON_INTEGRATION_POINTS_HELPER) #undef AKANTU_SPECIALIZE_INTEGRATE_ON_INTEGRATION_POINTS_HELPER #undef INTEGRATE } // namespace details } // namespace fe_engine template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>:: integrateOnIntegrationPoints(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 = getNbIntegrationPoints(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); fe_engine::details::IntegrateOnIntegrationPointsHelper<kind>::call( integrator, f, intf, nb_degree_of_freedom, type, ghost_type, filter_elements); } /* -------------------------------------------------------------------------- */ /** * Helper class to be able to write a partial specialization on the element kind */ namespace fe_engine { namespace details { template <ElementKind kind> struct InterpolateOnIntegrationPointsHelper { template <class S> static void call(const S &, const Array<Real> &, Array<Real> &, const UInt, const ElementType &, const GhostType &, const Array<UInt> &) { AKANTU_DEBUG_TO_IMPLEMENT(); } }; #define INTERPOLATE(type) \ shape_functions.template interpolateOnIntegrationPoints<type>( \ u, uq, nb_degree_of_freedom, ghost_type, filter_elements); #define AKANTU_SPECIALIZE_INTERPOLATE_ON_INTEGRATION_POINTS_HELPER(kind) \ template <> struct InterpolateOnIntegrationPointsHelper<kind> { \ template <class S> \ static void call(const S & shape_functions, const Array<Real> & u, \ Array<Real> & uq, const UInt nb_degree_of_freedom, \ const ElementType & type, const GhostType & ghost_type, \ const Array<UInt> & filter_elements) { \ AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTERPOLATE, kind); \ } \ }; AKANTU_BOOST_ALL_KIND_LIST( AKANTU_SPECIALIZE_INTERPOLATE_ON_INTEGRATION_POINTS_HELPER, AKANTU_FE_ENGINE_LIST_INTERPOLATE_ON_INTEGRATION_POINTS) #undef AKANTU_SPECIALIZE_INTERPOLATE_ON_INTEGRATION_POINTS_HELPER #undef INTERPOLATE } // namespace details } // namespace fe_engine template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>:: interpolateOnIntegrationPoints(const Array<Real> & u, Array<Real> & uq, const 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.getIntegrationPoints(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 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); fe_engine::details::InterpolateOnIntegrationPointsHelper<kind>::call( shape_functions, u, uq, nb_degree_of_freedom, type, ghost_type, filter_elements); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>:: interpolateOnIntegrationPoints( const Array<Real> & u, ElementTypeMapArray<Real> & uq, const ElementTypeMapArray<UInt> * filter_elements) const { AKANTU_DEBUG_IN(); const Array<UInt> * filter = NULL; for (auto ghost_type : ghost_types) { for (auto & type : uq.elementTypes(_all_dimensions, ghost_type, kind)) { UInt nb_quad_per_element = getNbIntegrationPoints(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); interpolateOnIntegrationPoints(u, quad, quad.getNbComponent(), type, ghost_type, *filter); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>:: computeIntegrationPointsCoordinates( ElementTypeMapArray<Real> & quadrature_points_coordinates, const ElementTypeMapArray<UInt> * filter_elements) const { const Array<Real> & nodes_coordinates = mesh.getNodes(); interpolateOnIntegrationPoints( nodes_coordinates, quadrature_points_coordinates, filter_elements); } /* -------------------------------------------------------------------------- */ template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>:: computeIntegrationPointsCoordinates( Array<Real> & quadrature_points_coordinates, const ElementType & type, const GhostType & ghost_type, const Array<UInt> & filter_elements) const { const Array<Real> & nodes_coordinates = mesh.getNodes(); UInt spatial_dimension = mesh.getSpatialDimension(); interpolateOnIntegrationPoints( nodes_coordinates, quadrature_points_coordinates, spatial_dimension, type, ghost_type, filter_elements); } /* -------------------------------------------------------------------------- */ template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>:: initElementalFieldInterpolationFromIntegrationPoints( const ElementTypeMapArray<Real> & interpolation_points_coordinates, ElementTypeMapArray<Real> & interpolation_points_coordinates_matrices, ElementTypeMapArray<Real> & quad_points_coordinates_inv_matrices, const ElementTypeMapArray<UInt> * element_filter) const { AKANTU_DEBUG_IN(); UInt spatial_dimension = this->mesh.getSpatialDimension(); ElementTypeMapArray<Real> quadrature_points_coordinates( "quadrature_points_coordinates_for_interpolation", getID(), getMemoryID()); quadrature_points_coordinates.initialize(*this, _nb_component = spatial_dimension); computeIntegrationPointsCoordinates(quadrature_points_coordinates, element_filter); shape_functions.initElementalFieldInterpolationFromIntegrationPoints( interpolation_points_coordinates, interpolation_points_coordinates_matrices, quad_points_coordinates_inv_matrices, quadrature_points_coordinates, element_filter); } /* -------------------------------------------------------------------------- */ template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>:: interpolateElementalFieldFromIntegrationPoints( const ElementTypeMapArray<Real> & field, const ElementTypeMapArray<Real> & interpolation_points_coordinates, ElementTypeMapArray<Real> & result, const GhostType ghost_type, const ElementTypeMapArray<UInt> * element_filter) const { ElementTypeMapArray<Real> interpolation_points_coordinates_matrices( "interpolation_points_coordinates_matrices", id, memory_id); ElementTypeMapArray<Real> quad_points_coordinates_inv_matrices( "quad_points_coordinates_inv_matrices", id, memory_id); initElementalFieldInterpolationFromIntegrationPoints( interpolation_points_coordinates, interpolation_points_coordinates_matrices, quad_points_coordinates_inv_matrices, element_filter); interpolateElementalFieldFromIntegrationPoints( field, interpolation_points_coordinates_matrices, quad_points_coordinates_inv_matrices, result, ghost_type, element_filter); } /* -------------------------------------------------------------------------- */ template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>:: interpolateElementalFieldFromIntegrationPoints( const ElementTypeMapArray<Real> & field, const ElementTypeMapArray<Real> & interpolation_points_coordinates_matrices, const ElementTypeMapArray<Real> & quad_points_coordinates_inv_matrices, ElementTypeMapArray<Real> & result, const GhostType ghost_type, const ElementTypeMapArray<UInt> * element_filter) const { shape_functions.interpolateElementalFieldFromIntegrationPoints( field, interpolation_points_coordinates_matrices, quad_points_coordinates_inv_matrices, result, ghost_type, element_filter); } /* -------------------------------------------------------------------------- */ /** * Helper class to be able to write a partial specialization on the element kind */ namespace fe_engine { namespace details { template <ElementKind kind> struct InterpolateHelper { template <class S> static void call(const S &, const Vector<Real> &, UInt, const Matrix<Real> &, Vector<Real> &, const ElementType &, const GhostType &) { AKANTU_DEBUG_TO_IMPLEMENT(); } }; #define INTERPOLATE(type) \ shape_functions.template interpolate<type>( \ real_coords, element, nodal_values, interpolated, ghost_type); #define AKANTU_SPECIALIZE_INTERPOLATE_HELPER(kind) \ template <> struct InterpolateHelper<kind> { \ template <class S> \ static void call(const S & shape_functions, \ const Vector<Real> & real_coords, UInt element, \ const Matrix<Real> & nodal_values, \ Vector<Real> & interpolated, const ElementType & type, \ const GhostType & ghost_type) { \ AKANTU_BOOST_KIND_ELEMENT_SWITCH(INTERPOLATE, kind); \ } \ }; AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_INTERPOLATE_HELPER, AKANTU_FE_ENGINE_LIST_INTERPOLATE) #undef AKANTU_SPECIALIZE_INTERPOLATE_HELPER #undef INTERPOLATE } // namespace details } // namespace fe_engine template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::interpolate( const Vector<Real> & real_coords, const Matrix<Real> & nodal_values, Vector<Real> & interpolated, const Element & element) const { AKANTU_DEBUG_IN(); fe_engine::details::InterpolateHelper<kind>::call( shape_functions, real_coords, element.element, nodal_values, interpolated, element.type, element.ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>:: computeNormalsOnIntegrationPoints(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); computeNormalsOnIntegrationPoints(mesh.getNodes(), ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>:: computeNormalsOnIntegrationPoints(const Array<Real> & field, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); // Real * coord = mesh.getNodes().storage(); UInt spatial_dimension = mesh.getSpatialDimension(); // allocate the normal arrays for (auto & type : mesh.elementTypes(element_dimension, ghost_type, kind)) { UInt size = mesh.getNbElement(type, ghost_type); if (normals_on_integration_points.exists(type, ghost_type)) { normals_on_integration_points(type, ghost_type).resize(size); } else { normals_on_integration_points.alloc(size, spatial_dimension, type, ghost_type); } } // loop over the type to build the normals for (auto & type : mesh.elementTypes(element_dimension, ghost_type, kind)) { Array<Real> & normals_on_quad = normals_on_integration_points(type, ghost_type); computeNormalsOnIntegrationPoints(field, normals_on_quad, type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Helper class to be able to write a partial specialization on the element kind */ namespace fe_engine { namespace details { template <ElementKind kind> struct ComputeNormalsOnIntegrationPoints { template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind k, class IOF> static void call(const FEEngineTemplate<I, S, k, IOF> &, const Array<Real> &, Array<Real> &, const ElementType &, const GhostType &) { AKANTU_DEBUG_TO_IMPLEMENT(); } }; #define COMPUTE_NORMALS_ON_INTEGRATION_POINTS(type) \ fem.template computeNormalsOnIntegrationPoints<type>(field, normal, \ ghost_type); #define AKANTU_SPECIALIZE_COMPUTE_NORMALS_ON_INTEGRATION_POINTS(kind) \ template <> struct ComputeNormalsOnIntegrationPoints<kind> { \ template <template <ElementKind, class> class I, \ template <ElementKind> class S, ElementKind k, class IOF> \ static void call(const FEEngineTemplate<I, S, k, IOF> & fem, \ const Array<Real> & field, Array<Real> & normal, \ const ElementType & type, const GhostType & ghost_type) { \ AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_NORMALS_ON_INTEGRATION_POINTS, \ kind); \ } \ }; AKANTU_BOOST_ALL_KIND_LIST( AKANTU_SPECIALIZE_COMPUTE_NORMALS_ON_INTEGRATION_POINTS, AKANTU_FE_ENGINE_LIST_COMPUTE_NORMALS_ON_INTEGRATION_POINTS) #undef AKANTU_SPECIALIZE_COMPUTE_NORMALS_ON_INTEGRATION_POINTS #undef COMPUTE_NORMALS_ON_INTEGRATION_POINTS } // namespace details } // namespace fe_engine template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>:: computeNormalsOnIntegrationPoints(const Array<Real> & field, Array<Real> & normal, const ElementType & type, const GhostType & ghost_type) const { fe_engine::details::ComputeNormalsOnIntegrationPoints<kind>::call( *this, field, normal, type, ghost_type); } /* -------------------------------------------------------------------------- */ template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> template <ElementType type> void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>:: computeNormalsOnIntegrationPoints(const Array<Real> & field, Array<Real> & normal, const GhostType & ghost_type) const { AKANTU_DEBUG_IN(); + if(type == _point_1) { + computeNormalsOnIntegrationPointsPoint1(field, normal, ghost_type); + return; + } + UInt spatial_dimension = mesh.getSpatialDimension(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_points = getNbIntegrationPoints(type, ghost_type); UInt nb_element = mesh.getConnectivity(type, ghost_type).getSize(); normal.resize(nb_element * nb_points); Array<Real>::matrix_iterator normals_on_quad = normal.begin_reinterpret(spatial_dimension, nb_points, nb_element); Array<Real> f_el(0, spatial_dimension * nb_nodes_per_element); FEEngine::extractNodalToElementField(mesh, field, f_el, type, ghost_type); const Matrix<Real> & quads = integrator.template getIntegrationPoints<type>(ghost_type); Array<Real>::matrix_iterator 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(); } /* -------------------------------------------------------------------------- */ /** * Helper class to be able to write a partial specialization on the element kind */ template <ElementKind kind> struct InverseMapHelper { template <class S> static void call(const S & shape_functions, const Vector<Real> & real_coords, UInt element, const ElementType & type, Vector<Real> & natural_coords, const GhostType & ghost_type) { AKANTU_DEBUG_TO_IMPLEMENT(); } }; #define INVERSE_MAP(type) \ shape_functions.template inverseMap<type>(real_coords, element, \ natural_coords, ghost_type); #define AKANTU_SPECIALIZE_INVERSE_MAP_HELPER(kind) \ template <> struct InverseMapHelper<kind> { \ template <class S> \ static void call(const S & shape_functions, \ const Vector<Real> & real_coords, UInt element, \ const ElementType & type, Vector<Real> & natural_coords, \ const GhostType & ghost_type) { \ AKANTU_BOOST_KIND_ELEMENT_SWITCH(INVERSE_MAP, kind); \ } \ }; AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_INVERSE_MAP_HELPER, AKANTU_FE_ENGINE_LIST_INVERSE_MAP) #undef AKANTU_SPECIALIZE_INVERSE_MAP_HELPER #undef INVERSE_MAP template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::inverseMap( const Vector<Real> & real_coords, UInt element, const ElementType & type, Vector<Real> & natural_coords, const GhostType & ghost_type) const { AKANTU_DEBUG_IN(); InverseMapHelper<kind>::call(shape_functions, real_coords, element, type, natural_coords, ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Helper class to be able to write a partial specialization on the element kind */ namespace fe_engine { namespace details { template <ElementKind kind> struct ContainsHelper { template <class S> static void call(const S &, const Vector<Real> &, UInt, const ElementType &, const GhostType &) { AKANTU_DEBUG_TO_IMPLEMENT(); } }; #define CONTAINS(type) \ contain = shape_functions.template contains<type>(real_coords, element, \ ghost_type); #define AKANTU_SPECIALIZE_CONTAINS_HELPER(kind) \ template <> struct ContainsHelper<kind> { \ template <template <ElementKind> class S, ElementKind k> \ static bool call(const S<k> & shape_functions, \ const Vector<Real> & real_coords, UInt element, \ const ElementType & type, const GhostType & ghost_type) { \ bool contain = false; \ AKANTU_BOOST_KIND_ELEMENT_SWITCH(CONTAINS, kind); \ return contain; \ } \ }; AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_CONTAINS_HELPER, AKANTU_FE_ENGINE_LIST_CONTAINS) #undef AKANTU_SPECIALIZE_CONTAINS_HELPER #undef CONTAINS } // namespace details } // namespace fe_engine template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> inline bool FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::contains( const Vector<Real> & real_coords, UInt element, const ElementType & type, const GhostType & ghost_type) const { return fe_engine::details::ContainsHelper<kind>::call( shape_functions, real_coords, element, type, ghost_type); } /* -------------------------------------------------------------------------- */ /** * Helper class to be able to write a partial specialization on the element kind */ namespace fe_engine { namespace details { template <ElementKind kind> struct ComputeShapesHelper { template <class S> static void call(const S &, const Vector<Real> &, UInt, const ElementType, Vector<Real> &, const GhostType &) { AKANTU_DEBUG_TO_IMPLEMENT(); } }; #define COMPUTE_SHAPES(type) \ shape_functions.template computeShapes<type>(real_coords, element, shapes, \ ghost_type); #define AKANTU_SPECIALIZE_COMPUTE_SHAPES_HELPER(kind) \ template <> struct ComputeShapesHelper<kind> { \ template <class S> \ static void call(const S & shape_functions, \ const Vector<Real> & real_coords, UInt element, \ const ElementType type, Vector<Real> & shapes, \ const GhostType & ghost_type) { \ AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_SHAPES, kind); \ } \ }; AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_COMPUTE_SHAPES_HELPER, AKANTU_FE_ENGINE_LIST_COMPUTE_SHAPES) #undef AKANTU_SPECIALIZE_COMPUTE_SHAPES_HELPER #undef COMPUTE_SHAPES } // namespace details } // namespace fe_engine template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::computeShapes( const Vector<Real> & real_coords, UInt element, const ElementType & type, Vector<Real> & shapes, const GhostType & ghost_type) const { AKANTU_DEBUG_IN(); fe_engine::details::ComputeShapesHelper<kind>::call( shape_functions, real_coords, element, type, shapes, ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Helper class to be able to write a partial specialization on the element kind */ namespace fe_engine { namespace details { template <ElementKind kind> struct ComputeShapeDerivativesHelper { template <class S> static void call(__attribute__((unused)) const S & shape_functions, __attribute__((unused)) const Vector<Real> & real_coords, __attribute__((unused)) UInt element, __attribute__((unused)) const ElementType type, __attribute__((unused)) Matrix<Real> & shape_derivatives, __attribute__((unused)) const GhostType & ghost_type) { AKANTU_DEBUG_TO_IMPLEMENT(); } }; #define COMPUTE_SHAPE_DERIVATIVES(type) \ Matrix<Real> coords_mat(real_coords.storage(), shape_derivatives.rows(), 1); \ Tensor3<Real> shapesd_tensor(shape_derivatives.storage(), \ shape_derivatives.rows(), \ shape_derivatives.cols(), 1); \ shape_functions.template computeShapeDerivatives<type>( \ coords_mat, element, shapesd_tensor, ghost_type); #define AKANTU_SPECIALIZE_COMPUTE_SHAPE_DERIVATIVES_HELPER(kind) \ template <> struct ComputeShapeDerivativesHelper<kind> { \ template <class S> \ static void call(const S & shape_functions, \ const Vector<Real> & real_coords, UInt element, \ const ElementType type, Matrix<Real> & shape_derivatives, \ const GhostType & ghost_type) { \ AKANTU_BOOST_KIND_ELEMENT_SWITCH(COMPUTE_SHAPE_DERIVATIVES, kind); \ } \ }; AKANTU_BOOST_ALL_KIND_LIST( AKANTU_SPECIALIZE_COMPUTE_SHAPE_DERIVATIVES_HELPER, AKANTU_FE_ENGINE_LIST_COMPUTE_SHAPES_DERIVATIVES) #undef AKANTU_SPECIALIZE_COMPUTE_SHAPE_DERIVATIVES_HELPER #undef COMPUTE_SHAPE_DERIVATIVES } // namespace details } // namespace fe_engine template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::computeShapeDerivatives( const Vector<Real> & real_coords, UInt element, const ElementType & type, Matrix<Real> & shape_derivatives, const GhostType & ghost_type) const { AKANTU_DEBUG_IN(); fe_engine::details::ComputeShapeDerivativesHelper<kind>::call( shape_functions, real_coords, element, type, shape_derivatives, ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Helper class to be able to write a partial specialization on the element kind */ namespace fe_engine { namespace details { template <ElementKind kind> struct GetNbIntegrationPointsHelper {}; #define GET_NB_INTEGRATION_POINTS(type) \ nb_quad_points = integrator.template getNbIntegrationPoints<type>(ghost_type); #define AKANTU_SPECIALIZE_GET_NB_INTEGRATION_POINTS_HELPER(kind) \ template <> struct GetNbIntegrationPointsHelper<kind> { \ template <template <ElementKind, class> class I, ElementKind k, class IOF> \ static UInt call(const I<k, IOF> & integrator, const ElementType type, \ const GhostType & ghost_type) { \ UInt nb_quad_points = 0; \ AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_NB_INTEGRATION_POINTS, kind); \ return nb_quad_points; \ } \ }; AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_GET_NB_INTEGRATION_POINTS_HELPER) #undef AKANTU_SPECIALIZE_GET_NB_INTEGRATION_POINTS_HELPER #undef GET_NB_INTEGRATION } // namespace details } // namespace fe_engine template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> inline UInt FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::getNbIntegrationPoints( const ElementType & type, const GhostType & ghost_type) const { return fe_engine::details::GetNbIntegrationPointsHelper<kind>::call( integrator, type, ghost_type); } /* -------------------------------------------------------------------------- */ /** * Helper class to be able to write a partial specialization on the element kind */ namespace fe_engine { namespace details { template <ElementKind kind> struct GetShapesHelper {}; #define GET_SHAPES(type) ret = &(shape_functions.getShapes(type, ghost_type)); #define AKANTU_SPECIALIZE_GET_SHAPES_HELPER(kind) \ template <> struct GetShapesHelper<kind> { \ template <class S> \ static const Array<Real> & call(const S & shape_functions, \ const ElementType type, \ const GhostType & ghost_type) { \ const Array<Real> * ret = NULL; \ AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_SHAPES, kind); \ return *ret; \ } \ }; AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_GET_SHAPES_HELPER) #undef AKANTU_SPECIALIZE_GET_SHAPES_HELPER #undef GET_SHAPES } // namespace details } // namespace fe_engine template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> inline const Array<Real> & FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::getShapes( const ElementType & type, const GhostType & ghost_type, __attribute__((unused)) UInt id) const { return fe_engine::details::GetShapesHelper<kind>::call(shape_functions, type, ghost_type); } /* -------------------------------------------------------------------------- */ /** * Helper class to be able to write a partial specialization on the element kind */ namespace fe_engine { namespace details { template <ElementKind kind> struct GetShapesDerivativesHelper { template <template <ElementKind> class S, ElementKind k> static const Array<Real> & call(const S<k> &, const ElementType &, const GhostType &, UInt) { AKANTU_DEBUG_TO_IMPLEMENT(); } }; #define GET_SHAPES_DERIVATIVES(type) \ ret = &(shape_functions.getShapesDerivatives(type, ghost_type)); #define AKANTU_SPECIALIZE_GET_SHAPES_DERIVATIVES_HELPER(kind) \ template <> struct GetShapesDerivativesHelper<kind> { \ template <template <ElementKind> class S, ElementKind k> \ static const Array<Real> & \ call(const S<k> & shape_functions, const ElementType type, \ const GhostType & ghost_type, __attribute__((unused)) UInt id) { \ const Array<Real> * ret = NULL; \ AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_SHAPES_DERIVATIVES, kind); \ return *ret; \ } \ }; AKANTU_BOOST_ALL_KIND_LIST(AKANTU_SPECIALIZE_GET_SHAPES_DERIVATIVES_HELPER, AKANTU_FE_ENGINE_LIST_GET_SHAPES_DERIVATIVES) #undef AKANTU_SPECIALIZE_GET_SHAPE_DERIVATIVES_HELPER #undef GET_SHAPES_DERIVATIVES } // namespace details } // namespace fe_engine template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> inline const Array<Real> & FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::getShapesDerivatives( const ElementType & type, const GhostType & ghost_type, __attribute__((unused)) UInt id) const { return fe_engine::details::GetShapesDerivativesHelper<kind>::call( shape_functions, type, ghost_type, id); } /* -------------------------------------------------------------------------- */ /** * Helper class to be able to write a partial specialization on the element kind */ namespace fe_engine { namespace details { template <ElementKind kind> struct GetIntegrationPointsHelper {}; #define GET_INTEGRATION_POINTS(type) \ ret = &(integrator.template getIntegrationPoints<type>(ghost_type)); #define AKANTU_SPECIALIZE_GET_INTEGRATION_POINTS_HELPER(kind) \ template <> struct GetIntegrationPointsHelper<kind> { \ template <template <ElementKind, class> class I, ElementKind k, class IOF> \ static const Matrix<Real> & call(const I<k, IOF> & integrator, \ const ElementType type, \ const GhostType & ghost_type) { \ const Matrix<Real> * ret = NULL; \ AKANTU_BOOST_KIND_ELEMENT_SWITCH(GET_INTEGRATION_POINTS, kind); \ return *ret; \ } \ }; AKANTU_BOOST_ALL_KIND(AKANTU_SPECIALIZE_GET_INTEGRATION_POINTS_HELPER) #undef AKANTU_SPECIALIZE_GET_INTEGRATION_POINTS_HELPER #undef GET_INTEGRATION_POINTS } // namespace details } // namespace fe_engine template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> inline const Matrix<Real> & FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::getIntegrationPoints( const ElementType & type, const GhostType & ghost_type) const { return fe_engine::details::GetIntegrationPointsHelper<kind>::call( integrator, type, ghost_type); } /* -------------------------------------------------------------------------- */ template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::printself( std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "FEEngineTemplate [" << std::endl; stream << space << " + parent [" << std::endl; FEEngine::printself(stream, indent + 3); stream << space << " ]" << std::endl; stream << space << " + shape functions [" << std::endl; shape_functions.printself(stream, indent + 3); stream << space << " ]" << std::endl; stream << space << " + integrator [" << std::endl; integrator.printself(stream, indent + 3); stream << space << " ]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::onElementsAdded( const Array<Element> & new_elements, const NewElementsEvent &) { integrator.onElementsAdded(new_elements); shape_functions.onElementsAdded(new_elements); } /* -------------------------------------------------------------------------- */ template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::onElementsRemoved( const Array<Element> &, const ElementTypeMapArray<UInt> &, const RemovedElementsEvent &) {} /* -------------------------------------------------------------------------- */ template <template <ElementKind, class> class I, template <ElementKind> class S, ElementKind kind, class IntegrationOrderFunctor> void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>::onElementsChanged( const Array<Element> &, const Array<Element> &, const ElementTypeMapArray<UInt> &, const ChangedElementsEvent &) {} /* -------------------------------------------------------------------------- */ - -} // namespace akantu - -#include "integrator_gauss.hh" -#include "shape_lagrange.hh" - -namespace akantu { - -/* -------------------------------------------------------------------------- */ -template <> -template <> -inline void FEEngineTemplate<IntegratorGauss, ShapeLagrange, _ek_regular, - DefaultIntegrationOrderFunctor>:: - computeNormalsOnIntegrationPoints<_point_1>( +template <template <ElementKind, class> class I, template <ElementKind> class S, + ElementKind kind, class IntegrationOrderFunctor> +inline void FEEngineTemplate<I, S, kind, IntegrationOrderFunctor>:: + computeNormalsOnIntegrationPointsPoint1( const Array<Real> &, 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 = getNbIntegrationPoints(type, ghost_type); UInt nb_element = mesh.getConnectivity(type, ghost_type).getSize(); normal.resize(nb_element * nb_points); Array<Real>::matrix_iterator normals_on_quad = normal.begin_reinterpret(spatial_dimension, nb_points, nb_element); const Array<std::vector<Element>> & segments = mesh.getElementToSubelement(type, ghost_type); const 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) { const 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(); } } // namespace akantu diff --git a/src/fe_engine/integrator_gauss.hh b/src/fe_engine/integrator_gauss.hh index 105a18f12..d35db8ed3 100644 --- a/src/fe_engine/integrator_gauss.hh +++ b/src/fe_engine/integrator_gauss.hh @@ -1,200 +1,199 @@ /** * @file integrator_gauss.hh * * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> * @author Nicolas Richart <nicolas.richart@epfl.ch> * * @date creation: Fri Jun 18 2010 * @date last modification: Thu Oct 22 2015 * * @brief Gauss integration facilities * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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_INTEGRATOR_GAUSS_HH__ -#define __AKANTU_INTEGRATOR_GAUSS_HH__ - /* -------------------------------------------------------------------------- */ #include "integrator.hh" /* -------------------------------------------------------------------------- */ +#ifndef __AKANTU_INTEGRATOR_GAUSS_HH__ +#define __AKANTU_INTEGRATOR_GAUSS_HH__ + namespace akantu { namespace integrator { namespace details { template <ElementKind> struct GaussIntegratorComputeJacobiansHelper; } // namespace details } // namespace fe_engine /* -------------------------------------------------------------------------- */ template <ElementKind kind, class IntegrationOrderFunctor> class IntegratorGauss : public Integrator { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: IntegratorGauss(const Mesh & mesh, const ID & id = "integrator_gauss", const MemoryID & memory_id = 0); virtual ~IntegratorGauss(){}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void initIntegrator(const Array<Real> & nodes, const ElementType & type, const GhostType & ghost_type); template <ElementType type> inline void initIntegrator(const Array<Real> & nodes, const GhostType & ghost_type); /// integrate f on the element "elem" of type "type" template <ElementType type> inline void integrateOnElement(const Array<Real> & f, Real * intf, UInt nb_degree_of_freedom, const UInt elem, const GhostType & ghost_type) const; /// integrate f for all elements of type "type" template <ElementType type> void integrate(const Array<Real> & in_f, Array<Real> & intf, UInt nb_degree_of_freedom, const GhostType & ghost_type, const Array<UInt> & filter_elements) const; /// integrate scalar field in_f template <ElementType type, UInt polynomial_degree> Real integrate(const Array<Real> & in_f, const GhostType & ghost_type = _not_ghost) const; /// integrate partially around a quadrature point (@f$ intf_q = f_q * J_q * /// w_q @f$) template <ElementType type> Real integrate(const Vector<Real> & in_f, UInt index, const GhostType & ghost_type) const; /// integrate scalar field in_f template <ElementType type> Real integrate(const Array<Real> & in_f, const GhostType & ghost_type, const Array<UInt> & filter_elements) const; /// integrate a field without using the pre-computed values template <ElementType type, UInt polynomial_degree> void integrate(const Array<Real> & in_f, Array<Real> & intf, UInt nb_degree_of_freedom, const GhostType & ghost_type) const; /// integrate partially around a quadrature point (@f$ intf_q = f_q * J_q * /// w_q @f$) template <ElementType type> void integrateOnIntegrationPoints(const Array<Real> & in_f, Array<Real> & intf, UInt nb_degree_of_freedom, const GhostType & ghost_type, const Array<UInt> & filter_elements) const; /// return a matrix with quadrature points natural coordinates template <ElementType type> const Matrix<Real> & getIntegrationPoints(const GhostType & ghost_type) const; /// return number of quadrature points template <ElementType type> UInt getNbIntegrationPoints(const GhostType & ghost_type) const; template <ElementType type, UInt n> Matrix<Real> getIntegrationPoints() const; template <ElementType type, UInt n> Vector<Real> getIntegrationWeights() const; protected: friend struct integrator::details::GaussIntegratorComputeJacobiansHelper<kind>; template <ElementType type> void computeJacobiansOnIntegrationPoints( const Array<Real> & nodes, const Matrix<Real> & quad_points, Array<Real> & jacobians, const GhostType & ghost_type, const Array<UInt> & filter_elements = empty_filter) const; void computeJacobiansOnIntegrationPoints( const Array<Real> & nodes, const Matrix<Real> & quad_points, Array<Real> & jacobians, const ElementType & type, const GhostType & ghost_type, const Array<UInt> & filter_elements = empty_filter) const; /// precompute jacobians on elements of type "type" template <ElementType type> void precomputeJacobiansOnQuadraturePoints(const Array<Real> & nodes, const GhostType & ghost_type); // multiply the jacobians by the integration weights and stores the results in // jacobians template <ElementType type, UInt polynomial_degree> void multiplyJacobiansByWeights(Array<Real> & jacobians) const; /// compute the vector of quadrature points natural coordinates template <ElementType type> void computeQuadraturePoints(const GhostType & ghost_type); /// check that the jacobians are not negative template <ElementType type> void checkJacobians(const GhostType & ghost_type) const; /// internal integrate partially around a quadrature point (@f$ intf_q = f_q * /// J_q * /// w_q @f$) void integrateOnIntegrationPoints(const Array<Real> & in_f, Array<Real> & intf, UInt nb_degree_of_freedom, const Array<Real> & jacobians, UInt nb_element) const; void integrate(const Array<Real> & in_f, Array<Real> & intf, UInt nb_degree_of_freedom, const Array<Real> & jacobians, UInt nb_element) const; public: /// compute the jacobians on quad points for a given element template <ElementType type> void computeJacobianOnQuadPointsByElement(const Matrix<Real> & node_coords, const Matrix<Real> & integration_points, Vector<Real> & jacobians) const; public: void onElementsAdded(const Array<Element> & elements) override; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// integrate the field f with the jacobian jac -> inte inline void integrate(Real * f, Real * jac, Real * inte, UInt nb_degree_of_freedom, UInt nb_quadrature_points) const; private: /// ElementTypeMap of the quadrature points ElementTypeMap<Matrix<Real>> quadrature_points; }; } // namespace akantu #include "integrator_gauss_inline_impl.cc" #endif /* __AKANTU_INTEGRATOR_GAUSS_HH__ */ diff --git a/src/fe_engine/shape_cohesive.hh b/src/fe_engine/shape_cohesive.hh index 4dc9d37f8..38e512944 100644 --- a/src/fe_engine/shape_cohesive.hh +++ b/src/fe_engine/shape_cohesive.hh @@ -1,170 +1,169 @@ /** * @file shape_cohesive.hh * * @author Nicolas Richart <nicolas.richart@epfl.ch> * @author Marco Vocialta <marco.vocialta@epfl.ch> * * @date creation: Tue Feb 15 2011 * @date last modification: Thu Oct 22 2015 * * @brief shape functions for cohesive elements description * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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_array.hh" -#include "shape_functions.hh" +#include "shape_lagrange.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SHAPE_COHESIVE_HH__ #define __AKANTU_SHAPE_COHESIVE_HH__ namespace akantu { struct CohesiveReduceFunctionMean { inline Real operator()(Real u_plus, Real u_minus) { return .5 * (u_plus + u_minus); } }; struct CohesiveReduceFunctionOpening { inline Real operator()(Real u_plus, Real u_minus) { return (u_plus - u_minus); } }; template <> class ShapeLagrange<_ek_cohesive> : public ShapeLagrangeBase { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: ShapeLagrange(const Mesh & mesh, const ID & id = "shape_cohesive", const MemoryID & memory_id = 0); virtual ~ShapeLagrange() = default; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: inline void initShapeFunctions(const Array<Real> & nodes, const Matrix<Real> & integration_points, const ElementType & type, const GhostType & ghost_type); /// extract the nodal values and store them per element template <ElementType type, class ReduceFunction> void extractNodalToElementField( const Array<Real> & nodal_f, Array<Real> & elemental_f, const GhostType & ghost_type = _not_ghost, const Array<UInt> & filter_elements = empty_filter) const; /// computes the shape functions derivatives for given interpolation points template <ElementType type> void computeShapeDerivativesOnIntegrationPoints( const Array<Real> & nodes, const Matrix<Real> & integration_points, Array<Real> & shape_derivatives, const GhostType & ghost_type, const Array<UInt> & filter_elements = empty_filter) const; void computeShapeDerivativesOnIntegrationPoints( const Array<Real> & nodes, const Matrix<Real> & integration_points, Array<Real> & shape_derivatives, const ElementType & type, const GhostType & ghost_type, const Array<UInt> & filter_elements) const override; /// pre compute all shapes on the element integration points from natural /// coordinates template <ElementType type> void precomputeShapesOnIntegrationPoints(const Array<Real> & nodes, GhostType ghost_type); /// pre compute all shape derivatives on the element integration points from /// natural coordinates template <ElementType type> void precomputeShapeDerivativesOnIntegrationPoints(const Array<Real> & nodes, GhostType ghost_type); /// interpolate nodal values on the integration points template <ElementType type, class ReduceFunction> void interpolateOnIntegrationPoints( const Array<Real> & u, Array<Real> & uq, UInt nb_degree_of_freedom, const GhostType ghost_type = _not_ghost, const Array<UInt> & filter_elements = empty_filter) const; template <ElementType type> void interpolateOnIntegrationPoints( const Array<Real> & u, Array<Real> & uq, UInt nb_degree_of_freedom, const GhostType ghost_type = _not_ghost, const Array<UInt> & filter_elements = empty_filter) const { interpolateOnIntegrationPoints<type, CohesiveReduceFunctionMean>( u, uq, nb_degree_of_freedom, ghost_type, filter_elements); } /// compute the gradient of u on the integration points in the natural /// coordinates template <ElementType type> void gradientOnIntegrationPoints( const Array<Real> & u, Array<Real> & nablauq, UInt nb_degree_of_freedom, GhostType ghost_type = _not_ghost, const Array<UInt> & filter_elements = empty_filter) const { variationOnIntegrationPoints<type, CohesiveReduceFunctionMean>( u, nablauq, nb_degree_of_freedom, ghost_type, filter_elements); } /// compute the gradient of u on the integration points template <ElementType type, class ReduceFunction> void variationOnIntegrationPoints( const Array<Real> & u, Array<Real> & nablauq, UInt nb_degree_of_freedom, GhostType ghost_type = _not_ghost, const Array<UInt> & filter_elements = empty_filter) const; /// compute the normals to the field u on integration points template <ElementType type, class ReduceFunction> void computeNormalsOnIntegrationPoints( const Array<Real> & u, Array<Real> & normals_u, GhostType ghost_type = _not_ghost, const Array<UInt> & filter_elements = empty_filter) const; /// multiply a field by shape functions template <ElementType type> void fieldTimesShapes(__attribute__((unused)) const Array<Real> & field, __attribute__((unused)) Array<Real> & fiedl_times_shapes, __attribute__((unused)) GhostType ghost_type) const { AKANTU_DEBUG_TO_IMPLEMENT(); } }; /// standard output stream operator template <class ShapeFunction> inline std::ostream & operator<<(std::ostream & stream, const ShapeCohesive<ShapeFunction> & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "shape_cohesive_inline_impl.cc" #endif /* __AKANTU_SHAPE_COHESIVE_HH__ */ diff --git a/src/fe_engine/shape_lagrange_base.cc b/src/fe_engine/shape_lagrange_base.cc index 29656b330..b857f8afb 100644 --- a/src/fe_engine/shape_lagrange_base.cc +++ b/src/fe_engine/shape_lagrange_base.cc @@ -1,123 +1,124 @@ /** * @file shape_lagrange_base.cc * * @author Nicolas Richart * * @date creation Thu Jul 27 2017 * * @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 "shape_lagrange_base.hh" +#include "mesh_iterators.hh" /* -------------------------------------------------------------------------- */ namespace akantu { ShapeLagrangeBase::ShapeLagrangeBase(const Mesh & mesh, const ElementKind & kind, const ID & id, const MemoryID & memory_id) : ShapeFunctions(mesh, id, memory_id), shapes("shapes_generic", id, memory_id), shapes_derivatives("shapes_derivatives_generic", id, memory_id), _kind(kind) {} /* -------------------------------------------------------------------------- */ ShapeLagrangeBase::~ShapeLagrangeBase() = default; /* -------------------------------------------------------------------------- */ void ShapeLagrangeBase::computeShapesOnIntegrationPoints( const Array<Real> & nodes, const Matrix<Real> & integration_points, Array<Real> & shapes, const ElementType & type, const GhostType & ghost_type, const Array<UInt> & filter_elements) const { #define AKANTU_COMPUTE_SHAPES(type) \ this->computeShapesOnIntegrationPoints<type>( \ nodes, integration_points, shapes, ghost_type, filter_elements) #define AKANTU_LAGRANGE_ELEMENT_TYPE \ AKANTU_ek_regular_ELEMENT_TYPE AKANTU_ek_cohesive_ELEMENT_TYPE AKANTU_BOOST_ELEMENT_SWITCH(AKANTU_COMPUTE_SHAPES, AKANTU_LAGRANGE_ELEMENT_TYPE); #undef AKANTU_COMPUTE_SHAPES //#undef AKANTU_LAGRANGE_ELEMENT_TYPE } /* -------------------------------------------------------------------------- */ void ShapeLagrangeBase::onElementsAdded(const Array<Element> & new_elements) { AKANTU_DEBUG_IN(); const auto & nodes = mesh.getNodes(); for (auto elements_range : MeshElementsByTypes(new_elements)) { auto type = elements_range.getType(); auto ghost_type = elements_range.getGhostType(); if (mesh.getKind(type) != _kind) continue; auto & elements = elements_range.getElements(); auto itp_type = FEEngine::getInterpolationType(type); if (not this->shapes_derivatives.exists(itp_type, ghost_type)) { auto size_of_shapesd = this->getShapeDerivativesSize(type); this->shapes_derivatives.alloc(0, size_of_shapesd, itp_type, ghost_type); } if (not shapes.exists(itp_type, ghost_type)) { auto size_of_shapes = this->getShapeSize(type); this->shapes.alloc(0, size_of_shapes, itp_type, ghost_type); } const auto & natural_coords = integration_points(type, ghost_type); computeShapesOnIntegrationPoints(nodes, natural_coords, shapes(itp_type, ghost_type), type, ghost_type, elements); computeShapeDerivativesOnIntegrationPoints( nodes, natural_coords, shapes_derivatives(itp_type, ghost_type), type, ghost_type, elements); } #undef INIT_SHAPE_FUNCTIONS AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ShapeLagrangeBase::onElementsRemoved( const Array<Element> &, const ElementTypeMapArray<UInt> & new_numbering) { this->shapes.onElementsRemoved(new_numbering); this->shapes_derivatives.onElementsRemoved(new_numbering); } /* -------------------------------------------------------------------------- */ void ShapeLagrangeBase::printself(std::ostream & stream, int indent) const { std::string space; for (Int i = 0; i < indent; i++, space += AKANTU_INDENT) ; stream << space << "Shapes Lagrange [" << std::endl; ShapeFunctions::printself(stream, indent + 1); shapes.printself(stream, indent + 1); shapes_derivatives.printself(stream, indent + 1); stream << space << "]" << std::endl; } } // namespace akantu diff --git a/src/mesh/element_group.hh b/src/mesh/element_group.hh index 181d9c3d6..69ff2de25 100644 --- a/src/mesh/element_group.hh +++ b/src/mesh/element_group.hh @@ -1,187 +1,188 @@ /** * @file element_group.hh * * @author Dana Christen <dana.christen@gmail.com> * @author Nicolas Richart <nicolas.richart@epfl.ch> * * @date creation: Fri May 03 2013 * @date last modification: Tue Aug 18 2015 * * @brief Stores information relevent to the notion of domain boundary and * surfaces. * * @section LICENSE * * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see <http://www.gnu.org/licenses/>. * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "dumpable.hh" #include "element_type_map.hh" #include "node_group.hh" /* -------------------------------------------------------------------------- */ #include <set> /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_ELEMENT_GROUP_HH__ #define __AKANTU_ELEMENT_GROUP_HH__ namespace akantu { class Mesh; class Element; } // namespace akantu namespace akantu { /* -------------------------------------------------------------------------- */ 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: using ElementList = ElementTypeMapArray<UInt>; using NodeList = Array<UInt>; /* ------------------------------------------------------------------------ */ /* Element iterator */ /* ------------------------------------------------------------------------ */ using type_iterator = ElementList::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; inline auto elementTypes(UInt dim = _all_dimensions, const GhostType & ghost_type = _not_ghost, const ElementKind & kind = _ek_regular) const { return elements.elementTypes(dim, ghost_type, kind); } using const_element_iterator = Array<UInt>::const_iterator<UInt>; inline const_element_iterator begin(const ElementType & type, const GhostType & ghost_type = _not_ghost) const; inline const_element_iterator end(const ElementType & type, 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); /// add an element to the group. By default the it does not add the nodes to /// the group inline void add(const Element & el, bool add_nodes = false, bool check_for_duplicate = true); /// \todo fix the default for add_nodes : make it coherent with the other /// method inline void add(const ElementType & type, UInt element, const GhostType & ghost_type = _not_ghost, bool add_nodes = true, bool check_for_duplicate = true); inline void addNode(UInt node_id, bool check_for_duplicate = true); inline void removeNode(UInt node_id); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; /// fill the elements based on the underlying node group. virtual void fillFromNodeGroup(); // sort and remove duplicated values void optimize(); private: inline void addElement(const ElementType & elem_type, UInt elem_id, const GhostType & ghost_type); friend class GroupManager; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: - AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Elements, elements, UInt); + const Array<UInt> & getElements(const ElementType & type, + const GhostType & ghost_type = _not_ghost) const; AKANTU_GET_MACRO(Elements, elements, const ElementTypeMapArray<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; /// empty arry for the iterator to work when an element type not present Array<UInt> empty_elements; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const ElementGroup & _this) { _this.printself(stream); return stream; } } // namespace akantu #include "element.hh" #include "element_group_inline_impl.cc" #endif /* __AKANTU_ELEMENT_GROUP_HH__ */ diff --git a/src/mesh/element_group_inline_impl.cc b/src/mesh/element_group_inline_impl.cc index 7ff809849..ebd1bd249 100644 --- a/src/mesh/element_group_inline_impl.cc +++ b/src/mesh/element_group_inline_impl.cc @@ -1,135 +1,146 @@ /** * @file element_group_inline_impl.cc * * @author Dana Christen <dana.christen@gmail.com> * @author Nicolas Richart <nicolas.richart@epfl.ch> * * @date creation: Wed Nov 13 2013 * @date last modification: Tue Aug 18 2015 * * @brief Stores information relevent to the notion of domain boundary and * surfaces. * * @section LICENSE * * Copyright (©) 2014, 2015 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 "mesh.hh" #include "element_group.hh" +#include "mesh.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_ELEMENT_GROUP_INLINE_IMPL_CC__ #define __AKANTU_ELEMENT_GROUP_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ inline void ElementGroup::add(const Element & el, bool add_nodes, bool check_for_duplicate) { this->add(el.type, el.element, el.ghost_type, add_nodes, check_for_duplicate); } /* -------------------------------------------------------------------------- */ inline void ElementGroup::add(const ElementType & type, UInt element, const GhostType & ghost_type, bool add_nodes, bool check_for_duplicate) { addElement(type, element, ghost_type); if (add_nodes) { Array<UInt>::const_vector_iterator it = mesh.getConnectivity(type, ghost_type) .begin(mesh.getNbNodesPerElement(type)) + element; const Vector<UInt> & conn = *it; for (UInt i = 0; i < conn.size(); ++i) addNode(conn[i], check_for_duplicate); } } /* -------------------------------------------------------------------------- */ inline void ElementGroup::addNode(UInt node_id, bool check_for_duplicate) { node_group.add(node_id, check_for_duplicate); } /* -------------------------------------------------------------------------- */ inline void ElementGroup::removeNode(UInt node_id) { node_group.remove(node_id); } /* -------------------------------------------------------------------------- */ inline void ElementGroup::addElement(const ElementType & elem_type, UInt elem_id, const GhostType & ghost_type) { if (!(elements.exists(elem_type, ghost_type))) { elements.alloc(0, 1, elem_type, ghost_type); } elements(elem_type, ghost_type).push_back(elem_id); this->dimension = UInt( std::max(Int(this->dimension), Int(mesh.getSpatialDimension(elem_type)))); } /* -------------------------------------------------------------------------- */ inline UInt ElementGroup::getNbNodes() const { return node_group.getSize(); } /* -------------------------------------------------------------------------- */ inline ElementGroup::type_iterator ElementGroup::firstType(UInt dim, const GhostType & ghost_type, const ElementKind & kind) const { return elements.firstType(dim, ghost_type, kind); } /* -------------------------------------------------------------------------- */ inline ElementGroup::type_iterator ElementGroup::lastType(UInt dim, const GhostType & ghost_type, const ElementKind & kind) const { return elements.lastType(dim, ghost_type, kind); } /* -------------------------------------------------------------------------- */ inline ElementGroup::const_element_iterator ElementGroup::begin(const ElementType & type, - const GhostType & ghost_type) const { + const GhostType & ghost_type) const { if (elements.exists(type, ghost_type)) { return elements(type, ghost_type).begin(); } else { return empty_elements.begin(); } } /* -------------------------------------------------------------------------- */ inline ElementGroup::const_element_iterator ElementGroup::end(const ElementType & type, - const GhostType & ghost_type) const { + const GhostType & ghost_type) const { if (elements.exists(type, ghost_type)) { return elements(type, ghost_type).end(); } else { return empty_elements.end(); } } +/* -------------------------------------------------------------------------- */ +inline const Array<UInt> & +ElementGroup::getElements(const ElementType & type, + const GhostType & ghost_type) const { + if (elements.exists(type, ghost_type)) { + return elements(type, ghost_type); + } else { + return empty_elements; + } +} + /* -------------------------------------------------------------------------- */ } // namespace akantu #endif /* __AKANTU_ELEMENT_GROUP_INLINE_IMPL_CC__ */ diff --git a/src/mesh/group_manager.hh b/src/mesh/group_manager.hh index 731f71f2c..caba3d89e 100644 --- a/src/mesh/group_manager.hh +++ b/src/mesh/group_manager.hh @@ -1,304 +1,305 @@ /** * @file group_manager.hh * * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> * @author Dana Christen <dana.christen@gmail.com> * @author David Simon Kammer <david.kammer@epfl.ch> * @author Nicolas Richart <nicolas.richart@epfl.ch> * @author Marco Vocialta <marco.vocialta@epfl.ch> * * @date creation: Wed Nov 13 2013 * @date last modification: Mon Nov 16 2015 * * @brief Stores information relevent to the notion of element and nodes *groups. * * @section LICENSE * * Copyright (©) 2014, 2015 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 "aka_common.hh" #include "element_type_map.hh" /* -------------------------------------------------------------------------- */ #include <set> /* -------------------------------------------------------------------------- */ namespace akantu { class ElementGroup; class NodeGroup; class Mesh; class Element; class ElementSynchronizer; template <bool> class CommunicationBufferTemplated; namespace dumper { class Field; } } // namespace akantu namespace akantu { /* -------------------------------------------------------------------------- */ class GroupManager { /* ------------------------------------------------------------------------ */ /* Typedefs */ /* ------------------------------------------------------------------------ */ private: using ElementGroups = std::map<std::string, ElementGroup *>; using NodeGroups = std::map<std::string, NodeGroup *>; public: typedef std::set<ElementType> GroupManagerTypeSet; /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: GroupManager(const Mesh & mesh, const ID & id = "group_manager", const MemoryID & memory_id = 0); virtual ~GroupManager(); /* ------------------------------------------------------------------------ */ /* Groups iterators */ /* ------------------------------------------------------------------------ */ 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; #ifndef SWIG #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); #endif /* ------------------------------------------------------------------------ */ /* Clustering filter */ /* ------------------------------------------------------------------------ */ public: class ClusteringFilter { public: virtual bool operator()(const Element &) const { return true; } }; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// create an empty node group NodeGroup & createNodeGroup(const std::string & group_name, bool replace_group = false); /// create a node group from another node group but filtered template <typename T> NodeGroup & createFilteredNodeGroup(const std::string & group_name, const NodeGroup & node_group, T & filter); /// 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 = _all_dimensions, bool replace_group = false); /// create an element group from another element group but filtered template <typename T> ElementGroup & createFilteredElementGroup(const std::string & group_name, UInt dimension, const NodeGroup & node_group, T & filter); /// destroy an element group and the associated node group void destroyElementGroup(const std::string & group_name, bool destroy_node_group = false); /// destroy all element groups and the associated node groups void destroyAllElementGroups(bool destroy_node_groups = false); /// create a element group using an existing node group ElementGroup & createElementGroup(const std::string & group_name, UInt dimension, NodeGroup & node_group); /// 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, Mesh & mesh_facets, std::string cluster_name_prefix = "cluster", const ClusteringFilter & filter = ClusteringFilter()); /// create element clusters for a given dimension UInt createClusters(UInt element_dimension, std::string cluster_name_prefix = "cluster", const ClusteringFilter & filter = ClusteringFilter()); private: /// create element clusters for a given dimension UInt createClusters(UInt element_dimension, std::string cluster_name_prefix, const ClusteringFilter & filter, Mesh & mesh_facets); public: /// Create an ElementGroup based on a NodeGroup void createElementGroupFromNodeGroup(const std::string & name, const std::string & node_group, UInt dimension = _all_dimensions); virtual void printself(std::ostream & stream, int indent = 0) const; /// this function insure that the group names are present on all processors /// /!\ it is a SMP call void synchronizeGroupNames(); /// register an elemental field to the given group name (overloading for /// ElementalPartionField) #ifndef SWIG template <typename T, template <bool> class dump_type> dumper::Field * createElementalField( const ElementTypeMapArray<T> & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap<UInt> nb_data_per_elem = ElementTypeMap<UInt>()); /// register an elemental field to the given group name (overloading for /// ElementalField) template <typename T, template <class> class ret_type, template <class, template <class> class, bool> class dump_type> dumper::Field * createElementalField( const ElementTypeMapArray<T> & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap<UInt> nb_data_per_elem = ElementTypeMap<UInt>()); /// register an elemental field to the given group name (overloading for /// MaterialInternalField) template <typename T, /// type of InternalMaterialField template <typename, bool filtered> class dump_type> dumper::Field * createElementalField(const ElementTypeMapArray<T> & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap<UInt> nb_data_per_elem); template <typename type, bool flag, template <class, bool> class ftype> dumper::Field * createNodalField(const ftype<type, flag> * field, const std::string & group_name, UInt padding_size = 0); template <typename type, bool flag, template <class, bool> class ftype> dumper::Field * createStridedNodalField(const ftype<type, flag> * field, const std::string & group_name, UInt size, UInt stride, UInt padding_size); protected: /// fill a buffer with all the group names void fillBufferWithGroupNames( CommunicationBufferTemplated<false> & comm_buffer) const; /// take a buffer and create the missing groups localy void checkAndAddGroups(CommunicationBufferTemplated<true> & buffer); /// register an elemental field to the given group name template <class dump_type, typename field_type> inline dumper::Field * createElementalField(const field_type & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap<UInt> nb_data_per_elem); /// register an elemental field to the given group name template <class dump_type, typename field_type> inline dumper::Field * createElementalFilteredField(const field_type & field, const std::string & group_name, UInt spatial_dimension, const ElementKind & kind, ElementTypeMap<UInt> nb_data_per_elem); #endif /* ------------------------------------------------------------------------ */ /* Accessor */ /* ------------------------------------------------------------------------ */ public: const ElementGroup & getElementGroup(const std::string & name) const; const NodeGroup & getNodeGroup(const std::string & name) const; ElementGroup & getElementGroup(const std::string & name); NodeGroup & getNodeGroup(const std::string & name); UInt getNbElementGroups(UInt dimension = _all_dimensions) const; + UInt getNbNodeGroups() { return node_groups.size(); }; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ 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 element 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; } } // namespace akantu #endif /* __AKANTU_GROUP_MANAGER_HH__ */ diff --git a/src/mesh/mesh_inline_impl.cc b/src/mesh/mesh_inline_impl.cc index 612d3fe5c..223265a07 100644 --- a/src/mesh/mesh_inline_impl.cc +++ b/src/mesh/mesh_inline_impl.cc @@ -1,681 +1,685 @@ /** * @file mesh_inline_impl.cc * * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> * @author Dana Christen <dana.christen@epfl.ch> * @author Nicolas Richart <nicolas.richart@epfl.ch> * @author Marco Vocialta <marco.vocialta@epfl.ch> * * @date creation: Thu Jul 15 2010 * @date last modification: Thu Jan 21 2016 * * @brief Implementation of the inline functions of the mesh class * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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 "mesh.hh" /* -------------------------------------------------------------------------- */ #if defined(AKANTU_COHESIVE_ELEMENT) #include "cohesive_element.hh" #endif #ifndef __AKANTU_MESH_INLINE_IMPL_CC__ #define __AKANTU_MESH_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ inline RemovedNodesEvent::RemovedNodesEvent(const Mesh & mesh) : new_numbering(mesh.getNbNodes(), 1, "new_numbering") {} /* -------------------------------------------------------------------------- */ inline RemovedElementsEvent::RemovedElementsEvent(const Mesh & mesh, ID new_numbering_id) : new_numbering(new_numbering_id, mesh.getID(), mesh.getMemoryID()) {} /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent<NewElementsEvent>(NewElementsEvent & event) { this->nodes_to_elements.resize(nodes->getSize()); for (const auto & elem : event.getList()) { const Array<UInt> & conn = connectivities(elem.type, elem.ghost_type); UInt nb_nodes_per_elem = this->getNbNodesPerElement(elem.type); for (UInt n = 0; n < nb_nodes_per_elem; ++n) { UInt node = conn(elem.element, n); if (not nodes_to_elements[node]) nodes_to_elements[node] = std::make_unique<std::set<Element>>(); nodes_to_elements[node]->insert(elem); } } EventHandlerManager<MeshEventHandler>::sendEvent(event); } /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent<NewNodesEvent>(NewNodesEvent & event) { this->computeBoundingBox(); EventHandlerManager<MeshEventHandler>::sendEvent(event); } /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent<RemovedElementsEvent>(RemovedElementsEvent & event) { this->connectivities.onElementsRemoved(event.getNewNumbering()); this->fillNodesToElements(); this->computeBoundingBox(); EventHandlerManager<MeshEventHandler>::sendEvent(event); } /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent<RemovedNodesEvent>(RemovedNodesEvent & event) { const auto & new_numbering = event.getNewNumbering(); this->removeNodesFromArray(*nodes, new_numbering); if (nodes_global_ids) this->removeNodesFromArray(*nodes_global_ids, new_numbering); if (nodes_type.getSize() != 0) this->removeNodesFromArray(nodes_type, new_numbering); - std::vector<std::unique_ptr<std::set<Element>>> tmp(nodes_to_elements.size()); - auto it = nodes_to_elements.begin(); - - UInt new_nb_nodes = 0; - for (auto new_i : new_numbering) { - if (new_i != UInt(-1)) { - tmp[new_i] = std::move(*it); - ++new_nb_nodes; + if(not nodes_to_elements.empty()) { + std::vector<std::unique_ptr<std::set<Element>>> tmp(nodes_to_elements.size()); + auto it = nodes_to_elements.begin(); + + UInt new_nb_nodes = 0; + for (auto new_i : new_numbering) { + if (new_i != UInt(-1)) { + tmp[new_i] = std::move(*it); + ++new_nb_nodes; + } + ++it; } - ++it; - } - tmp.resize(new_nb_nodes); - std::move(tmp.begin(), tmp.end(), nodes_to_elements.begin()); + tmp.resize(new_nb_nodes); + std::move(tmp.begin(), tmp.end(), nodes_to_elements.begin()); + } computeBoundingBox(); EventHandlerManager<MeshEventHandler>::sendEvent(event); } /* -------------------------------------------------------------------------- */ template <typename T> inline void Mesh::removeNodesFromArray(Array<T> & vect, const Array<UInt> & new_numbering) { Array<T> tmp(vect.getSize(), vect.getNbComponent()); UInt nb_component = vect.getNbComponent(); UInt new_nb_nodes = 0; for (UInt i = 0; i < new_numbering.getSize(); ++i) { UInt new_i = new_numbering(i); if (new_i != UInt(-1)) { T * to_copy = vect.storage() + i * nb_component; std::uninitialized_copy(to_copy, to_copy + nb_component, tmp.storage() + new_i * nb_component); ++new_nb_nodes; } } tmp.resize(new_nb_nodes); vect.copy(tmp); } /* -------------------------------------------------------------------------- */ // inline UInt Mesh::elementToLinearized(const Element & elem) const { // AKANTU_DEBUG_ASSERT(elem.type < _max_element_type && // elem.element < types_offsets.storage()[elem.type + // 1], // "The element " << elem << "does not exists in the mesh // " // << getID()); // return types_offsets.storage()[elem.type] + elem.element; // } // /* -------------------------------------------------------------------------- // */ inline Element Mesh::linearizedToElement(UInt linearized_element) const { // UInt t; // for (t = _not_defined; // t != _max_element_type && linearized_element >= types_offsets(t); ++t) // ; // AKANTU_DEBUG_ASSERT(linearized_element < types_offsets(t), // "The linearized element " // << linearized_element // << "does not exists in the mesh " << getID()); // --t; // ElementType type = ElementType(t); // return Element(type, linearized_element - types_offsets.storage()[t], // _not_ghost, getKind(type)); // } // /* -------------------------------------------------------------------------- // */ inline void Mesh::updateTypesOffsets(const GhostType & ghost_type) { // Array<UInt> * types_offsets_ptr = &this->types_offsets; // if (ghost_type == _ghost) // types_offsets_ptr = &this->ghost_types_offsets; // Array<UInt> & types_offsets = *types_offsets_ptr; // types_offsets.clear(); // for (auto type : elementTypes(_all_dimensions, ghost_type, // _ek_not_defined)) // types_offsets(type) = connectivities(type, ghost_type).getSize(); // for (UInt t = _not_defined + 1; t < _max_element_type; ++t) // types_offsets(t) += types_offsets(t - 1); // for (UInt t = _max_element_type; t > _not_defined; --t) // types_offsets(t) = types_offsets(t - 1); // types_offsets(0) = 0; // } // /* -------------------------------------------------------------------------- // */ inline const Mesh::ConnectivityTypeList & // Mesh::getConnectivityTypeList(const GhostType & ghost_type) const { // if (ghost_type == _not_ghost) // return type_set; // else // return ghost_type_set; // } /* -------------------------------------------------------------------------- */ inline Array<UInt> & Mesh::getNodesGlobalIdsPointer() { AKANTU_DEBUG_IN(); if (not nodes_global_ids) { auto nb_nodes = nodes->getSize(); nodes_global_ids = std::make_unique<Array<UInt>>( nb_nodes, 1, getID() + ":nodes_global_ids"); UInt i = 0; for (auto & global_id : *nodes_global_ids) { global_id = ++i; } } AKANTU_DEBUG_OUT(); return *nodes_global_ids; } /* -------------------------------------------------------------------------- */ inline Array<NodeType> & Mesh::getNodesTypePointer() { AKANTU_DEBUG_IN(); if (nodes_type.getSize() == 0) { nodes_type.resize(nodes->getSize()); nodes_type.set(_nt_normal); } AKANTU_DEBUG_OUT(); return nodes_type; } /* -------------------------------------------------------------------------- */ inline Array<UInt> & Mesh::getConnectivityPointer(const ElementType & type, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); Array<UInt> * tmp; if (!connectivities.exists(type, ghost_type)) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); tmp = &(connectivities.alloc(0, nb_nodes_per_element, type, ghost_type)); AKANTU_DEBUG_INFO("The connectivity vector for the type " << type << " created"); // if (ghost_type == _not_ghost) // type_set.insert(type); // else // ghost_type_set.insert(type); // updateTypesOffsets(ghost_type); } else { tmp = &connectivities(type, ghost_type); } AKANTU_DEBUG_OUT(); return *tmp; } /* -------------------------------------------------------------------------- */ inline Array<std::vector<Element>> & Mesh::getElementToSubelementPointer(const ElementType & type, const GhostType & ghost_type) { return getDataPointer<std::vector<Element>>("element_to_subelement", type, ghost_type, 1, true); } /* -------------------------------------------------------------------------- */ inline Array<Element> & Mesh::getSubelementToElementPointer(const ElementType & type, const GhostType & ghost_type) { return getDataPointer<Element>("subelement_to_element", type, ghost_type, getNbFacetsPerElement(type), true, is_mesh_facets); } /* -------------------------------------------------------------------------- */ inline const Array<std::vector<Element>> & Mesh::getElementToSubelement(const ElementType & type, const GhostType & ghost_type) const { return getData<std::vector<Element>>("element_to_subelement", type, ghost_type); } /* -------------------------------------------------------------------------- */ inline Array<std::vector<Element>> & Mesh::getElementToSubelement(const ElementType & type, const GhostType & ghost_type) { return getData<std::vector<Element>>("element_to_subelement", type, ghost_type); } /* -------------------------------------------------------------------------- */ inline const Array<Element> & Mesh::getSubelementToElement(const ElementType & type, const GhostType & ghost_type) const { return getData<Element>("subelement_to_element", type, ghost_type); } /* -------------------------------------------------------------------------- */ inline Array<Element> & Mesh::getSubelementToElement(const ElementType & type, const GhostType & ghost_type) { return getData<Element>("subelement_to_element", type, ghost_type); } /* -------------------------------------------------------------------------- */ template <typename T> inline Array<T> & Mesh::getDataPointer(const std::string & data_name, const ElementType & el_type, const GhostType & ghost_type, UInt nb_component, bool size_to_nb_element, bool resize_with_parent) { Array<T> & tmp = mesh_data.getElementalDataArrayAlloc<T>( data_name, el_type, ghost_type, nb_component); if (size_to_nb_element) { if (resize_with_parent) tmp.resize(mesh_parent->getNbElement(el_type, ghost_type)); else tmp.resize(this->getNbElement(el_type, ghost_type)); } else { tmp.resize(0); } return tmp; } /* -------------------------------------------------------------------------- */ template <typename T> inline const Array<T> & Mesh::getData(const std::string & data_name, const ElementType & el_type, const GhostType & ghost_type) const { return mesh_data.getElementalDataArray<T>(data_name, el_type, ghost_type); } /* -------------------------------------------------------------------------- */ template <typename T> inline Array<T> & Mesh::getData(const std::string & data_name, const ElementType & el_type, const GhostType & ghost_type) { return mesh_data.getElementalDataArray<T>(data_name, el_type, ghost_type); } /* -------------------------------------------------------------------------- */ template <typename T> inline const ElementTypeMapArray<T> & Mesh::getData(const std::string & data_name) const { return mesh_data.getElementalData<T>(data_name); } /* -------------------------------------------------------------------------- */ template <typename T> inline ElementTypeMapArray<T> & Mesh::getData(const std::string & data_name) { return mesh_data.getElementalData<T>(data_name); } /* -------------------------------------------------------------------------- */ template <typename T> inline ElementTypeMapArray<T> & Mesh::registerData(const std::string & data_name) { this->mesh_data.registerElementalData<T>(data_name); return this->getData<T>(data_name); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbElement(const ElementType & type, const GhostType & ghost_type) const { try { const Array<UInt> & conn = connectivities(type, ghost_type); return conn.getSize(); } catch (...) { return 0; } } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbElement(const UInt spatial_dimension, const GhostType & ghost_type, const ElementKind & kind) const { AKANTU_DEBUG_ASSERT(spatial_dimension <= 3 || spatial_dimension == UInt(-1), "spatial_dimension is " << spatial_dimension << " and is greater than 3 !"); UInt nb_element = 0; for (auto type : elementTypes(spatial_dimension, ghost_type, kind)) nb_element += getNbElement(type, ghost_type); return nb_element; } /* -------------------------------------------------------------------------- */ inline void Mesh::getBarycenter(UInt element, const ElementType & type, Real * barycenter, GhostType ghost_type) const { UInt * conn_val = getConnectivity(type, ghost_type).storage(); UInt nb_nodes_per_element = getNbNodesPerElement(type); Real * local_coord = new Real[spatial_dimension * nb_nodes_per_element]; UInt offset = element * nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) { memcpy(local_coord + n * spatial_dimension, nodes->storage() + conn_val[offset + n] * spatial_dimension, spatial_dimension * sizeof(Real)); } Math::barycenter(local_coord, nb_nodes_per_element, spatial_dimension, barycenter); delete[] local_coord; } /* -------------------------------------------------------------------------- */ inline void Mesh::getBarycenter(const Element & element, Vector<Real> & barycenter) const { getBarycenter(element.element, element.type, barycenter.storage(), element.ghost_type); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbNodesPerElement(const ElementType & type) { UInt nb_nodes_per_element = 0; #define GET_NB_NODES_PER_ELEMENT(type) \ nb_nodes_per_element = ElementClass<type>::getNbNodesPerElement() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_NODES_PER_ELEMENT); #undef GET_NB_NODES_PER_ELEMENT return nb_nodes_per_element; } /* -------------------------------------------------------------------------- */ inline ElementType Mesh::getP1ElementType(const ElementType & type) { ElementType p1_type = _not_defined; #define GET_P1_TYPE(type) p1_type = ElementClass<type>::getP1ElementType() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_P1_TYPE); #undef GET_P1_TYPE return p1_type; } /* -------------------------------------------------------------------------- */ inline ElementKind Mesh::getKind(const ElementType & type) { ElementKind kind = _ek_not_defined; #define GET_KIND(type) kind = ElementClass<type>::getKind() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_KIND); #undef GET_KIND return kind; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getSpatialDimension(const ElementType & type) { UInt spatial_dimension = 0; #define GET_SPATIAL_DIMENSION(type) \ spatial_dimension = ElementClass<type>::getSpatialDimension() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_SPATIAL_DIMENSION); #undef GET_SPATIAL_DIMENSION return spatial_dimension; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbFacetTypes(const ElementType & type, __attribute__((unused)) UInt t) { UInt nb = 0; #define GET_NB_FACET_TYPE(type) nb = ElementClass<type>::getNbFacetTypes() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_FACET_TYPE); #undef GET_NB_FACET_TYPE return nb; } /* -------------------------------------------------------------------------- */ inline ElementType Mesh::getFacetType(const ElementType & type, UInt t) { ElementType surface_type = _not_defined; #define GET_FACET_TYPE(type) surface_type = ElementClass<type>::getFacetType(t) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_FACET_TYPE); #undef GET_FACET_TYPE return surface_type; } /* -------------------------------------------------------------------------- */ inline VectorProxy<ElementType> Mesh::getAllFacetTypes(const ElementType & type) { #define GET_FACET_TYPE(type) \ UInt nb = ElementClass<type>::getNbFacetTypes(); \ ElementType * elt_ptr = \ const_cast<ElementType *>(ElementClass<type>::getFacetTypeInternal()); \ return VectorProxy<ElementType>(elt_ptr, nb); AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_FACET_TYPE); #undef GET_FACET_TYPE } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbFacetsPerElement(const ElementType & type) { AKANTU_DEBUG_IN(); UInt n_facet = 0; #define GET_NB_FACET(type) n_facet = ElementClass<type>::getNbFacetsPerElement() AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_FACET); #undef GET_NB_FACET AKANTU_DEBUG_OUT(); return n_facet; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbFacetsPerElement(const ElementType & type, UInt t) { AKANTU_DEBUG_IN(); UInt n_facet = 0; #define GET_NB_FACET(type) \ n_facet = ElementClass<type>::getNbFacetsPerElement(t) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_NB_FACET); #undef GET_NB_FACET AKANTU_DEBUG_OUT(); return n_facet; } /* -------------------------------------------------------------------------- */ inline MatrixProxy<UInt> Mesh::getFacetLocalConnectivity(const ElementType & type, UInt t) { AKANTU_DEBUG_IN(); #define GET_FACET_CON(type) \ AKANTU_DEBUG_OUT(); \ return ElementClass<type>::getFacetLocalConnectivityPerElement(t) AKANTU_BOOST_ALL_ELEMENT_SWITCH(GET_FACET_CON); #undef GET_FACET_CON AKANTU_DEBUG_OUT(); return MatrixProxy<UInt>( nullptr, 0, 0); // This avoid a compilation warning but will certainly // also cause a segfault if reached } /* -------------------------------------------------------------------------- */ inline Matrix<UInt> Mesh::getFacetConnectivity(const Element & element, UInt t) const { AKANTU_DEBUG_IN(); Matrix<UInt> local_facets(getFacetLocalConnectivity(element.type, t), false); Matrix<UInt> facets(local_facets.rows(), local_facets.cols()); const Array<UInt> & conn = connectivities(element.type, element.ghost_type); for (UInt f = 0; f < facets.rows(); ++f) { for (UInt n = 0; n < facets.cols(); ++n) { facets(f, n) = conn(element.element, local_facets(f, n)); } } AKANTU_DEBUG_OUT(); return facets; } /* -------------------------------------------------------------------------- */ template <typename T> inline void Mesh::extractNodalValuesFromElement( const Array<T> & nodal_values, T * local_coord, UInt * connectivity, UInt n_nodes, UInt nb_degree_of_freedom) const { for (UInt n = 0; n < n_nodes; ++n) { memcpy(local_coord + n * nb_degree_of_freedom, nodal_values.storage() + connectivity[n] * nb_degree_of_freedom, nb_degree_of_freedom * sizeof(T)); } } /* -------------------------------------------------------------------------- */ inline void Mesh::addConnectivityType(const ElementType & type, const GhostType & ghost_type) { getConnectivityPointer(type, ghost_type); } /* -------------------------------------------------------------------------- */ inline bool Mesh::isPureGhostNode(UInt n) const { return nodes_type.getSize() ? (nodes_type(n) == _nt_pure_gost) : false; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isLocalOrMasterNode(UInt n) const { return nodes_type.getSize() ? (nodes_type(n) == _nt_master) || (nodes_type(n) == _nt_normal) : true; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isLocalNode(UInt n) const { return nodes_type.getSize() ? nodes_type(n) == _nt_normal : true; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isMasterNode(UInt n) const { return nodes_type.getSize() ? nodes_type(n) == _nt_master : false; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isSlaveNode(UInt n) const { return nodes_type.getSize() ? nodes_type(n) >= 0 : false; } /* -------------------------------------------------------------------------- */ inline NodeType Mesh::getNodeType(UInt local_id) const { return nodes_type.getSize() ? nodes_type(local_id) : _nt_normal; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNodeGlobalId(UInt local_id) const { return nodes_global_ids ? (*nodes_global_ids)(local_id) : local_id; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNodeLocalId(UInt global_id) const { AKANTU_DEBUG_ASSERT(nodes_global_ids != NULL, "The global ids are note set."); return nodes_global_ids->find(global_id); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbGlobalNodes() const { return nodes_global_ids ? nb_global_nodes : nodes->getSize(); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbNodesPerElementList(const Array<Element> & elements) { UInt nb_nodes_per_element = 0; UInt nb_nodes = 0; ElementType current_element_type = _not_defined; Array<Element>::const_iterator<Element> el_it = elements.begin(); Array<Element>::const_iterator<Element> el_end = elements.end(); for (; el_it != el_end; ++el_it) { const Element & el = *el_it; if (el.type != current_element_type) { current_element_type = el.type; nb_nodes_per_element = Mesh::getNbNodesPerElement(current_element_type); } nb_nodes += nb_nodes_per_element; } return nb_nodes; } /* -------------------------------------------------------------------------- */ inline Mesh & Mesh::getMeshFacets() { if (!this->mesh_facets) AKANTU_EXCEPTION( "No facet mesh is defined yet! check the buildFacets functions"); return *this->mesh_facets; } /* -------------------------------------------------------------------------- */ inline const Mesh & Mesh::getMeshFacets() const { if (!this->mesh_facets) AKANTU_EXCEPTION( "No facet mesh is defined yet! check the buildFacets functions"); return *this->mesh_facets; } /* -------------------------------------------------------------------------- */ inline const Mesh & Mesh::getMeshParent() const { if (!this->mesh_parent) AKANTU_EXCEPTION( "No parent mesh is defined! This is only valid in a mesh_facets"); return *this->mesh_parent; } /* -------------------------------------------------------------------------- */ } // namespace akantu #endif /* __AKANTU_MESH_INLINE_IMPL_CC__ */ diff --git a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.hh b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.hh index 981ec2cd5..69a3cae41 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.hh +++ b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.hh @@ -1,266 +1,265 @@ /** * @file material_cohesive.hh * * @author Nicolas Richart <nicolas.richart@epfl.ch> * @author Seyedeh Mohadeseh Taheri Mousavi <mohadeseh.taherimousavi@epfl.ch> * @author Marco Vocialta <marco.vocialta@epfl.ch> * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Jan 12 2016 * * @brief Specialization of the material class for cohesive elements * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see <http://www.gnu.org/licenses/>. * */ /* -------------------------------------------------------------------------- */ - -#include "aka_common.hh" #include "material.hh" -#include "fe_engine_template.hh" -#include "aka_common.hh" -#include "cohesive_internal_field.hh" +/* -------------------------------------------------------------------------- */ #include "cohesive_element_inserter.hh" - +#include "cohesive_internal_field.hh" +#include "fe_engine_template.hh" +#include "integrator_gauss.hh" +#include "shape_lagrange.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_COHESIVE_HH__ #define __AKANTU_MATERIAL_COHESIVE_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class SolidMechanicsModelCohesive; } namespace akantu { class MaterialCohesive : public Material { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef FEEngineTemplate<IntegratorGauss, ShapeLagrange, _ek_cohesive> MyFEEngineCohesiveType; public: MaterialCohesive(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialCohesive(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize the material computed parameter virtual void initMaterial(); /// compute tractions (including normals and openings) void computeTraction(GhostType ghost_type = _not_ghost); /// assemble residual void assembleInternalForces(GhostType ghost_type = _not_ghost); /// compute reversible and total energies by element void computeEnergies(); /// check stress for cohesive elements' insertion, by default it /// also updates the cohesive elements' data virtual void checkInsertion(__attribute__((unused)) bool check_only = false) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// check delta_max for cohesive elements in case of no convergence /// in the solveStep (only for extrinsic-implicit) virtual void checkDeltaMax(__attribute__((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// reset variables when convergence is reached (only for /// extrinsic-implicit) virtual void resetVariables(__attribute__((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// interpolate stress on given positions for each element (empty /// implemantation to avoid the generic call to be done on cohesive elements) virtual void interpolateStress(__attribute__((unused)) const ElementType type, __attribute__((unused)) Array<Real> & result){}; /// compute the stresses virtual void computeAllStresses(__attribute__((unused)) GhostType ghost_type = _not_ghost){}; // add the facet to be handled by the material UInt addFacet(const Element & element); protected: virtual void computeTangentTraction(__attribute__((unused)) const ElementType & el_type, __attribute__((unused)) Array<Real> & tangent_matrix, __attribute__((unused)) const Array<Real> & normal, __attribute__((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute the normal void computeNormal(const Array<Real> & position, Array<Real> & normal, ElementType type, GhostType ghost_type); /// compute the opening void computeOpening(const Array<Real> & displacement, Array<Real> & opening, ElementType type, GhostType ghost_type); template <ElementType type> void computeNormal(const Array<Real> & position, Array<Real> & normal, GhostType ghost_type); /// assemble stiffness void assembleStiffnessMatrix(GhostType ghost_type); /// constitutive law virtual void computeTraction(const Array<Real> & normal, ElementType el_type, GhostType ghost_type = _not_ghost) = 0; /// parallelism functions inline UInt getNbDataForElements(const Array<Element> & elements, SynchronizationTag tag) const; inline void packElementData(CommunicationBuffer & buffer, const Array<Element> & elements, SynchronizationTag tag) const; inline void unpackElementData(CommunicationBuffer & buffer, const Array<Element> & elements, SynchronizationTag tag); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the opening AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Opening, opening, Real); /// get the traction AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Traction, tractions, Real); /// get damage AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Damage, damage, Real); /// get facet filter AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(FacetFilter, facet_filter, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(FacetFilter, facet_filter, UInt); AKANTU_GET_MACRO(FacetFilter, facet_filter, const ElementTypeMapArray<UInt> &); // AKANTU_GET_MACRO(ElementFilter, element_filter, const // ElementTypeMapArray<UInt> &); /// compute reversible energy Real getReversibleEnergy(); /// compute dissipated energy Real getDissipatedEnergy(); /// compute contact energy Real getContactEnergy(); /// get energy virtual Real getEnergy(std::string type); /// return the energy (identified by id) for the provided element virtual Real getEnergy(std::string energy_id, ElementType type, UInt index) { return Material::getEnergy(energy_id, type, index); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// list of facets assigned to this material ElementTypeMapArray<UInt> facet_filter; /// Link to the cohesive fem object in the model MyFEEngineCohesiveType * fem_cohesive; private: /// reversible energy by quadrature point CohesiveInternalField<Real> reversible_energy; /// total energy by quadrature point CohesiveInternalField<Real> total_energy; protected: /// opening in all elements and quadrature points CohesiveInternalField<Real> opening; /// opening in all elements and quadrature points (previous time step) CohesiveInternalField<Real> opening_old; /// traction in all elements and quadrature points CohesiveInternalField<Real> tractions; /// traction in all elements and quadrature points (previous time step) CohesiveInternalField<Real> tractions_old; /// traction due to contact CohesiveInternalField<Real> contact_tractions; /// normal openings for contact tractions CohesiveInternalField<Real> contact_opening; /// maximum displacement CohesiveInternalField<Real> delta_max; /// tell if the previous delta_max state is needed (in iterative schemes) bool use_previous_delta_max; /// tell if the previous opening state is needed (in iterative schemes) bool use_previous_opening; /// damage CohesiveInternalField<Real> damage; /// pointer to the solid mechanics model for cohesive elements SolidMechanicsModelCohesive * model; /// critical stress RandomInternalField<Real, FacetInternalField> sigma_c; /// critical displacement Real delta_c; /// array to temporarily store the normals Array<Real> normal; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_cohesive_inline_impl.cc" -} // akantu +} // namespace akantu #include "cohesive_internal_field_tmpl.hh" #endif /* __AKANTU_MATERIAL_COHESIVE_HH__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc index f083dc51a..9570d37a6 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model.cc @@ -1,1261 +1,1267 @@ /** * @file solid_mechanics_model.cc * * @author Ramin Aghababaei <ramin.aghababaei@epfl.ch> * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> * @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch> * @author David Simon Kammer <david.kammer@epfl.ch> * @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch> * @author Nicolas Richart <nicolas.richart@epfl.ch> * @author Clement Roux <clement.roux@epfl.ch> * @author Marco Vocialta <marco.vocialta@epfl.ch> * * @date creation: Tue Jul 27 2010 * @date last modification: Tue Jan 19 2016 * * @brief Implementation of the SolidMechanicsModel class * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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 "solid_mechanics_model.hh" +#include "integrator_gauss.hh" +#include "shape_lagrange.hh" #include "element_synchronizer.hh" #include "sparse_matrix.hh" #include "static_communicator.hh" #include "synchronizer_registry.hh" #include "dumpable_inline_impl.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_paraview.hh" #endif #include "material_non_local.hh" /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ /** * A solid mechanics model need a mesh and a dimension to be created. the model * by it self can not do a lot, the good init functions should be called in * order to configure the model depending on what we want to do. * * @param mesh mesh representing the model we want to simulate * @param dim spatial dimension of the problem, if dim = 0 (default value) the * dimension of the problem is assumed to be the on of the mesh * @param id an id to identify the model */ SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : Model(mesh, dim, id, memory_id), BoundaryCondition<SolidMechanicsModel>(), f_m2a(1.0), displacement(nullptr), previous_displacement(nullptr), displacement_increment(nullptr), mass(nullptr), velocity(nullptr), acceleration(nullptr), external_force(nullptr), internal_force(nullptr), blocked_dofs(nullptr), current_position(nullptr), mass_matrix(nullptr), velocity_damping_matrix(nullptr), stiffness_matrix(nullptr), jacobian_matrix(nullptr), material_index("material index", id, memory_id), material_local_numbering("material local numbering", id, memory_id), material_selector(new DefaultMaterialSelector(material_index)), is_default_material_selector(true), increment_flag(false), are_materials_instantiated(false) { //, pbc_synch(nullptr) { AKANTU_DEBUG_IN(); this->registerFEEngineObject<MyFEEngineType>("SolidMechanicsFEEngine", mesh, Model::spatial_dimension); this->mesh.registerEventHandler(*this, _ehp_solid_mechanics_model); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper<DumperParaview>("paraview_all", id, true); this->mesh.addDumpMesh(mesh, Model::spatial_dimension, _not_ghost, _ek_regular); #endif this->initDOFManager(); this->registerDataAccessor(*this); if (this->mesh.isDistributed()) { auto & synchronizer = this->mesh.getElementSynchronizer(); this->registerSynchronizer(synchronizer, _gst_material_id); this->registerSynchronizer(synchronizer, _gst_smm_mass); this->registerSynchronizer(synchronizer, _gst_smm_stress); this->registerSynchronizer(synchronizer, _gst_smm_boundary); this->registerSynchronizer(synchronizer, _gst_for_dump); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModel::~SolidMechanicsModel() { AKANTU_DEBUG_IN(); if (is_default_material_selector) { delete material_selector; material_selector = nullptr; } for (auto & internal : this->registered_internals) { delete internal.second; } // delete pbc_synch; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setTimeStep(Real time_step, const ID & solver_id) { Model::setTimeStep(time_step, solver_id); #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper().setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ /* Initialization */ /* -------------------------------------------------------------------------- */ /** * This function groups many of the initialization in on function. For most of * basics case the function should be enough. The functions initialize the * model, the internal vectors, set them to 0, and depending on the parameters * it also initialize the explicit or implicit solver. * * @param material_file the file containing the materials to use * @param method the analysis method wanted. See the akantu::AnalysisMethod for * the different possibilities */ void SolidMechanicsModel::initFull(const ModelOptions & options) { Model::initFull(options); const SolidMechanicsModelOptions & smm_options = dynamic_cast<const SolidMechanicsModelOptions &>(options); this->method = smm_options.analysis_method; // initialize the vectors this->initArrays(); if (!this->hasDefaultSolver()) this->initNewSolver(this->method); // initialize pbc if (this->pbc_pair.size() != 0) this->initPBC(); // initialize the materials if (this->parser->getLastParsedFile() != "") { this->instantiateMaterials(); } //if (!smm_options.no_init_materials) { this->initMaterials(); //} // if (increment_flag) this->initBC(*this, *displacement, *displacement_increment, *external_force); // else // this->initBC(*this, *displacement, *external_force); } /* -------------------------------------------------------------------------- */ TimeStepSolverType SolidMechanicsModel::getDefaultSolverType() const { return _tsst_dynamic_lumped; } /* -------------------------------------------------------------------------- */ ModelSolverOptions SolidMechanicsModel::getDefaultSolverOptions( const TimeStepSolverType & type) const { ModelSolverOptions options; switch (type) { case _tsst_dynamic_lumped: { options.non_linear_solver_type = _nls_lumped; options.integration_scheme_type["displacement"] = _ist_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; break; } case _tsst_static: { options.non_linear_solver_type = _nls_newton_raphson; options.integration_scheme_type["displacement"] = _ist_pseudo_time; options.solution_type["displacement"] = IntegrationScheme::_not_defined; break; } case _tsst_dynamic: { if (this->method == _explicit_consistent_mass) { options.non_linear_solver_type = _nls_newton_raphson; options.integration_scheme_type["displacement"] = _ist_central_difference; options.solution_type["displacement"] = IntegrationScheme::_acceleration; } else { options.non_linear_solver_type = _nls_newton_raphson; options.integration_scheme_type["displacement"] = _ist_trapezoidal_rule_2; options.solution_type["displacement"] = IntegrationScheme::_displacement; } break; } } return options; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initNewSolver(const AnalysisMethod & method) { ID solver_name; TimeStepSolverType tss_type; this->method = method; switch (this->method) { case _explicit_lumped_mass: { solver_name = "explicit_lumped"; tss_type = _tsst_dynamic_lumped; break; } case _explicit_consistent_mass: { solver_name = "explicit"; tss_type = _tsst_dynamic; break; } case _static: { solver_name = "static"; tss_type = _tsst_static; break; } case _implicit_dynamic: { solver_name = "implicit"; tss_type = _tsst_dynamic; break; } } if (!this->hasSolver(solver_name)) { ModelSolverOptions options = this->getDefaultSolverOptions(tss_type); this->getNewSolver(solver_name, tss_type, options.non_linear_solver_type); this->setIntegrationScheme(solver_name, "displacement", options.integration_scheme_type["displacement"], options.solution_type["displacement"]); this->setDefaultSolver(solver_name); } } /* -------------------------------------------------------------------------- */ template <typename T> void SolidMechanicsModel::allocNodalField(Array<T> *& array, const ID & name) { if (array == nullptr) { UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_disp; sstr_disp << id << ":" << name; array = &(alloc<T>(sstr_disp.str(), nb_nodes, Model::spatial_dimension, T())); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initSolver( TimeStepSolverType time_step_solver_type, __attribute__((unused)) NonLinearSolverType non_linear_solver_type) { DOFManager & dof_manager = this->getDOFManager(); /* ------------------------------------------------------------------------ */ // for alloc type of solvers this->allocNodalField(this->displacement, "displacement"); this->allocNodalField(this->previous_displacement, "previous_displacement"); this->allocNodalField(this->displacement_increment, "displacement_increment"); this->allocNodalField(this->internal_force, "internal_force"); this->allocNodalField(this->external_force, "external_force"); this->allocNodalField(this->blocked_dofs, "blocked_dofs"); this->allocNodalField(this->current_position, "current_position"); // initialize the current positions this->current_position->copy(this->mesh.getNodes()); /* ------------------------------------------------------------------------ */ if (!dof_manager.hasDOFs("displacement")) { dof_manager.registerDOFs("displacement", *this->displacement, _dst_nodal); dof_manager.registerBlockedDOFs("displacement", *this->blocked_dofs); dof_manager.registerDOFsIncrement("displacement", *this->displacement_increment); dof_manager.registerDOFsPrevious("displacement", *this->previous_displacement); } /* ------------------------------------------------------------------------ */ // for dynamic if (time_step_solver_type == _tsst_dynamic || time_step_solver_type == _tsst_dynamic_lumped) { this->allocNodalField(this->velocity, "velocity"); this->allocNodalField(this->acceleration, "acceleration"); if (!dof_manager.hasDOFsDerivatives("displacement", 1)) { dof_manager.registerDOFsDerivative("displacement", 1, *this->velocity); dof_manager.registerDOFsDerivative("displacement", 2, *this->acceleration); } } if (time_step_solver_type == _tsst_dynamic || time_step_solver_type == _tsst_static) { if (!dof_manager.hasMatrix("K")) { dof_manager.getNewMatrix("K", _symmetric); } if (!dof_manager.hasMatrix("J")) { dof_manager.getNewMatrix("J", "K"); } } } /* -------------------------------------------------------------------------- */ // void SolidMechanicsModel::initParallel(MeshPartition & partition, // DataAccessor<Element> * data_accessor) // { // AKANTU_DEBUG_IN(); // if (data_accessor == nullptr) // 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); // synch_registry->registerSynchronizer(*synch_parallel, _gst_for_dump); // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initFEEngineBoundary() { FEEngine & fem_boundary = getFEEngineBoundary(); fem_boundary.initShapeFunctions(_not_ghost); fem_boundary.initShapeFunctions(_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_ghost); } /* -------------------------------------------------------------------------- */ /** * Allocate all the needed vectors. By default their are not necessarily set to * 0 */ void SolidMechanicsModel::initArrays() { AKANTU_DEBUG_IN(); for (auto ghost_type : ghost_types) { for (auto type : mesh.elementTypes(Model::spatial_dimension, ghost_type, _ek_not_defined)) { UInt nb_element = mesh.getNbElement(type, ghost_type); this->material_index.alloc(nb_element, 1, type, ghost_type); this->material_local_numbering.alloc(nb_element, 1, type, ghost_type); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian * */ void SolidMechanicsModel::initModel() { /// \todo add the current position as a parameter to initShapeFunctions for /// large deformation getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ // void SolidMechanicsModel::initPBC() { // Model::initPBC(); // registerPBCSynchronizer(); // // as long as there are ones on the diagonal of the matrix, we can put // // boudandary true for slaves // std::map<UInt, UInt>::iterator it = pbc_pair.begin(); // std::map<UInt, UInt>::iterator end = pbc_pair.end(); // UInt dim = mesh.getSpatialDimension(); // while (it != end) { // for (UInt i = 0; i < dim; ++i) // (*blocked_dofs)((*it).first, i) = true; // ++it; // } // } // /* -------------------------------------------------------------------------- // */ // void SolidMechanicsModel::registerPBCSynchronizer() { // pbc_synch = new PBCSynchronizer(pbc_pair); // synch_registry->registerSynchronizer(*pbc_synch, _gst_smm_uv); // synch_registry->registerSynchronizer(*pbc_synch, _gst_smm_mass); // synch_registry->registerSynchronizer(*pbc_synch, _gst_smm_res); // synch_registry->registerSynchronizer(*pbc_synch, _gst_for_dump); // // changeLocalEquationNumberForPBC(pbc_pair, mesh.getSpatialDimension()); // } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleResidual() { AKANTU_DEBUG_IN(); /* ------------------------------------------------------------------------ */ // computes the internal forces this->assembleInternalForces(); /* ------------------------------------------------------------------------ */ this->getDOFManager().assembleToResidual("displacement", *this->external_force, 1); this->getDOFManager().assembleToResidual("displacement", *this->internal_force, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleJacobian() { this->assembleStiffnessMatrix(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::predictor() { ++displacement_release; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::corrector() { ++displacement_release; } /* -------------------------------------------------------------------------- */ /** * This function computes the internal forces as F_{int} = \int_{\Omega} N * \sigma d\Omega@f$ */ void SolidMechanicsModel::assembleInternalForces() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); this->internal_force->clear(); // compute the stresses of local elements AKANTU_DEBUG_INFO("Compute local stresses"); for (auto & material : materials) { material->computeAllStresses(_not_ghost); } #ifdef AKANTU_DAMAGE_NON_LOCAL /* ------------------------------------------------------------------------ */ /* Computation of the non local part */ if(this->non_local_manager) this->non_local_manager->computeAllNonLocalStresses(); #endif // communicate the stresses AKANTU_DEBUG_INFO("Send data for residual assembly"); this->asynchronousSynchronize(_gst_smm_stress); // assemble the forces due to local stresses AKANTU_DEBUG_INFO("Assemble residual for local elements"); for (auto & material : materials) { material->assembleInternalForces(_not_ghost); } // finalize communications AKANTU_DEBUG_INFO("Wait distant stresses"); this->waitEndSynchronize(_gst_smm_stress); // assemble the stresses due to ghost elements AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); for (auto & material : materials) { material->assembleInternalForces(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix."); // Check if materials need to recompute the matrix bool need_to_reassemble = false; for (auto & material : materials) { need_to_reassemble |= material->hasStiffnessMatrixChanged(); } if (need_to_reassemble) { this->getDOFManager().getMatrix("K").clear(); // call compute stiffness matrix on each local elements for (auto & material : materials) { material->assembleStiffnessMatrix(_not_ghost); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateCurrentPosition() { if (this->current_position_release == this->displacement_release) return; this->current_position->copy(this->mesh.getNodes()); auto cpos_it = this->current_position->begin(Model::spatial_dimension); auto cpos_end = this->current_position->end(Model::spatial_dimension); auto disp_it = this->displacement->begin(Model::spatial_dimension); for (; cpos_it != cpos_end; ++cpos_it, ++disp_it) { *cpos_it += *disp_it; } this->current_position_release = this->displacement_release; } /* -------------------------------------------------------------------------- */ const Array<Real> & SolidMechanicsModel::getCurrentPosition() { this->updateCurrentPosition(); return *this->current_position; } /* -------------------------------------------------------------------------- */ // void SolidMechanicsModel::initializeUpdateResidualData() { // AKANTU_DEBUG_IN(); // UInt nb_nodes = mesh.getNbNodes(); // internal_force->resize(nb_nodes); // /// copy the forces in residual for boundary conditions // this->getDOFManager().assembleToResidual("displacement", // *this->external_force); // // start synchronization // this->asynchronousSynchronize(_gst_smm_uv); // this->waitEndSynchronize(_gst_smm_uv); // this->updateCurrentPosition(); // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ /* Explicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ // void SolidMechanicsModel::computeStresses() { // if (isExplicit()) { // // start synchronization // this->asynchronousSynchronize(_gst_smm_uv); // this->waitEndSynchronize(_gst_smm_uv); // // compute stresses on all local elements for each materials // std::vector<Material *>::iterator mat_it; // for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { // Material & mat = **mat_it; // mat.computeAllStresses(_not_ghost); // } // #ifdef AKANTU_DAMAGE_NON_LOCAL // /* Computation of the non local part */ // this->non_local_manager->computeAllNonLocalStresses(); // #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); // } // } // } /* -------------------------------------------------------------------------- */ /* Implicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ // /** // * Initialize the solver and create the sparse matrices needed. // * // */ // void SolidMechanicsModel::initSolver(__attribute__((unused)) // SolverOptions & options) { // UInt nb_global_nodes = mesh.getNbGlobalNodes(); // jacobian_matrix = &(this->getDOFManager().getNewMatrix("jacobian", // _symmetric)); // // jacobian_matrix->buildProfile(mesh, *dof_synchronizer, // spatial_dimension); // if (!isExplicit()) { // delete stiffness_matrix; // std::stringstream sstr_sti; // sstr_sti << id << ":stiffness_matrix"; // stiffness_matrix = &(this->getDOFManager().getNewMatrix("stiffness", // _symmetric)); // } // if (solver) solver->initialize(options); // } // /* -------------------------------------------------------------------------- // */ // void SolidMechanicsModel::initJacobianMatrix() { // // @todo make it more flexible: this is an ugly patch to treat the case of // non // // fix profile of the K matrix // delete jacobian_matrix; // jacobian_matrix = &(this->getDOFManager().getNewMatrix("jacobian", // "stiffness")); // std::stringstream sstr_solv; sstr_solv << id << ":solver"; // delete solver; // solver = new SolverMumps(*jacobian_matrix, sstr_solv.str()); // if(solver) // solver->initialize(_solver_no_options); // } /* -------------------------------------------------------------------------- */ /** * Initialize the implicit solver, either for dynamic or static cases, * * @param dynamic */ // void SolidMechanicsModel::initImplicit(bool dynamic) { // AKANTU_DEBUG_IN(); // method = dynamic ? _implicit_dynamic : _static; // if (!increment) // setIncrementFlagOn(); // initSolver(); // // if(method == _implicit_dynamic) { // // if(integrator) delete integrator; // // integrator = new TrapezoidalRule2(); // // } // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ // SparseMatrix & SolidMechanicsModel::initVelocityDampingMatrix() { // return this->getDOFManager().getNewMatrix("velocity_damping", "jacobian"); // } // /* -------------------------------------------------------------------------- // */ // void SolidMechanicsModel::implicitPred() { // AKANTU_DEBUG_IN(); // if(previous_displacement) previous_displacement->copy(*displacement); // if(method == _implicit_dynamic) // integrator->integrationSchemePred(time_step, // *displacement, // *velocity, // *acceleration, // *blocked_dofs); // AKANTU_DEBUG_OUT(); // } // /* -------------------------------------------------------------------------- // */ // void SolidMechanicsModel::implicitCorr() { // AKANTU_DEBUG_IN(); // if(method == _implicit_dynamic) { // integrator->integrationSchemeCorrDispl(time_step, // *displacement, // *velocity, // *acceleration, // *blocked_dofs, // *increment); // } else { // UInt nb_nodes = displacement->getSize(); // UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes; // Real * incr_val = increment->storage(); // Real * disp_val = displacement->storage(); // bool * boun_val = blocked_dofs->storage(); // for (UInt j = 0; j < nb_degree_of_freedom; ++j, ++disp_val, ++incr_val, // ++boun_val){ // *incr_val *= (1. - *boun_val); // *disp_val += *incr_val; // } // } // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ // void SolidMechanicsModel::updateIncrement() { // AKANTU_DEBUG_IN(); // auto incr_it = this->displacement_increment->begin(spatial_dimension); // auto incr_end = this->displacement_increment->end(spatial_dimension); // auto disp_it = this->displacement->begin(spatial_dimension); // auto prev_disp_it = this->previous_displacement->begin(spatial_dimension); // for (; incr_it != incr_end; ++incr_it) { // *incr_it = *disp_it; // *incr_it -= *prev_disp_it; // } // 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(); // } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateDataForNonLocalCriterion( ElementTypeMapReal & criterion) { const ID field_name = criterion.getName(); for (auto & material : materials) { if (!material->isInternal<Real>(field_name, _ek_regular)) continue; for (auto ghost_type : ghost_types) { material->flattenInternal(field_name, criterion, ghost_type, _ek_regular); } } } /* -------------------------------------------------------------------------- */ /* Information */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ // void SolidMechanicsModel::synchronizeBoundaries() { // AKANTU_DEBUG_IN(); // this->synchronize(_gst_smm_boundary); // AKANTU_DEBUG_OUT(); // } // /* -------------------------------------------------------------------------- // */ void SolidMechanicsModel::synchronizeResidual() { // AKANTU_DEBUG_IN(); // this->synchronize(_gst_smm_res); // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ // void SolidMechanicsModel::setIncrementFlagOn() { // AKANTU_DEBUG_IN(); // if (!displacement_increment) { // this->allocNodalField(displacement_increment, spatial_dimension, // "displacement_increment"); // this->getDOFManager().registerDOFsIncrement("displacement", // *this->displacement_increment); // } // 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, _so_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); Real min_dt = std::numeric_limits<Real>::max(); this->updateCurrentPosition(); Element elem; elem.ghost_type = ghost_type; elem.kind = _ek_regular; for (auto type : mesh.elementTypes(Model::spatial_dimension, ghost_type, _ek_regular)) { elem.type = type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); UInt nb_element = mesh.getNbElement(type); auto mat_indexes = material_index(type, ghost_type).begin(); auto mat_loc_num = material_local_numbering(type, ghost_type).begin(); Array<Real> X(0, nb_nodes_per_element * Model::spatial_dimension); FEEngine::extractNodalToElementField(mesh, *current_position, X, type, _not_ghost); auto X_el = X.begin(Model::spatial_dimension, nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++X_el, ++mat_indexes, ++mat_loc_num) { elem.element = *mat_loc_num; Real el_h = getFEEngine().getElementInradius(*X_el, type); Real el_c = this->materials[*mat_indexes]->getCelerity(elem); Real el_dt = el_h / el_c; min_dt = std::min(min_dt, el_dt); } } AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy() { AKANTU_DEBUG_IN(); Real ekin = 0.; UInt nb_nodes = mesh.getNbNodes(); if (this->getDOFManager().hasLumpedMatrix("M")) { auto m_it = this->mass->begin(Model::spatial_dimension); auto m_end = this->mass->end(Model::spatial_dimension); auto v_it = this->velocity->begin(Model::spatial_dimension); for (UInt n = 0; m_it != m_end; ++n, ++m_it, ++v_it) { const Vector<Real> & v = *v_it; const Vector<Real> & m = *m_it; Real mv2 = 0; bool is_local_node = mesh.isLocalOrMasterNode(n); // bool is_not_pbc_slave_node = !isPBCSlaveNode(n); bool count_node = is_local_node; // && is_not_pbc_slave_node; if (count_node) { for (UInt i = 0; i < Model::spatial_dimension; ++i) { if (m(i) > std::numeric_limits<Real>::epsilon()) mv2 += v(i) * v(i) * m(i); } } ekin += mv2; } } else if (this->getDOFManager().hasMatrix("M")) { Array<Real> Mv(nb_nodes, Model::spatial_dimension); this->getDOFManager().getMatrix("M").matVecMul(*this->velocity, Mv); auto mv_it = Mv.begin(Model::spatial_dimension); auto mv_end = Mv.end(Model::spatial_dimension); auto v_it = this->velocity->begin(Model::spatial_dimension); for (; mv_it != mv_end; ++mv_it, ++v_it) { ekin += v_it->dot(*mv_it); } } else { AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix."); } StaticCommunicator::getStaticCommunicator().allReduce(ekin, _so_sum); AKANTU_DEBUG_OUT(); return ekin * .5; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy(const ElementType & type, UInt index) { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type); Array<Real> vel_on_quad(nb_quadrature_points, Model::spatial_dimension); Array<UInt> filter_element(1, 1, index); getFEEngine().interpolateOnIntegrationPoints(*velocity, vel_on_quad, Model::spatial_dimension, type, _not_ghost, filter_element); Array<Real>::vector_iterator vit = vel_on_quad.begin(Model::spatial_dimension); Array<Real>::vector_iterator vend = vel_on_quad.end(Model::spatial_dimension); Vector<Real> rho_v2(nb_quadrature_points); Real rho = materials[material_index(type)(index)]->getRho(); for (UInt q = 0; vit != vend; ++vit, ++q) { rho_v2(q) = rho * vit->dot(*vit); } AKANTU_DEBUG_OUT(); return .5 * getFEEngine().integrate(rho_v2, type, index); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getExternalWork() { AKANTU_DEBUG_IN(); auto incr_or_velo_it = this->velocity->begin(Model::spatial_dimension); if (this->method == _static) { incr_or_velo_it = this->displacement_increment->begin(Model::spatial_dimension); } auto ext_force_it = external_force->begin(Model::spatial_dimension); auto int_force_it = internal_force->begin(Model::spatial_dimension); auto boun_it = blocked_dofs->begin(Model::spatial_dimension); Real work = 0.; UInt nb_nodes = this->mesh.getNbNodes(); for (UInt n = 0; n < nb_nodes; ++n, ++ext_force_it, ++int_force_it, ++boun_it, ++incr_or_velo_it) { const auto & int_force = *int_force_it; const auto & ext_force = *ext_force_it; const auto & boun = *boun_it; const auto & incr_or_velo = *incr_or_velo_it; bool is_local_node = this->mesh.isLocalOrMasterNode(n); // bool is_not_pbc_slave_node = !this->isPBCSlaveNode(n); bool count_node = is_local_node; // && is_not_pbc_slave_node; if (count_node) { for (UInt i = 0; i < Model::spatial_dimension; ++i) { if (boun(i)) work -= int_force(i) * incr_or_velo(i); else work += ext_force(i) * incr_or_velo(i); } } } StaticCommunicator::getStaticCommunicator().allReduce(work, _so_sum); if (this->method != _static) work *= this->getTimeStep(); 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.; for (auto & material : materials) energy += material->getEnergy(energy_id); /// reduction sum over all processors StaticCommunicator::getStaticCommunicator().allReduce(energy, _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); } UInt mat_index = this->material_index(type, _not_ghost)(index); UInt mat_loc_num = this->material_local_numbering(type, _not_ghost)(index); Real energy = this->materials[mat_index]->getEnergy(energy_id, type, mat_loc_num); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsAdded(const Array<Element> & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); for (auto ghost_type : ghost_types) { for (auto type : mesh.elementTypes(Model::spatial_dimension, ghost_type, _ek_not_defined)) { UInt nb_element = this->mesh.getNbElement(type, ghost_type); if (!material_index.exists(type, ghost_type)) { this->material_index.alloc(nb_element, 1, type, ghost_type); this->material_local_numbering.alloc(nb_element, 1, type, ghost_type); } else { this->material_index(type, ghost_type).resize(nb_element); this->material_local_numbering(type, ghost_type).resize(nb_element); } } } ElementTypeMapArray<UInt> filter("new_element_filter", this->getID(), this->getMemoryID()); for (auto & elem : element_list) { if (!filter.exists(elem.type, elem.ghost_type)) filter.alloc(0, 1, elem.type, elem.ghost_type); filter(elem.type, elem.ghost_type).push_back(elem.element); } this->assignMaterialToElements(&filter); for (auto & material : materials) material->onElementsAdded(element_list, event); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsRemoved( __attribute__((unused)) const Array<Element> & element_list, const ElementTypeMapArray<UInt> & new_numbering, const RemovedElementsEvent & event) { this->getFEEngine().initShapeFunctions(_not_ghost); this->getFEEngine().initShapeFunctions(_ghost); for (auto & material : materials) { material->onElementsRemoved(element_list, new_numbering, event); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesAdded(const Array<UInt> & nodes_list, __attribute__((unused)) const NewNodesEvent & event) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); if (displacement) { displacement->resize(nb_nodes, 0.); ++displacement_release; } if (mass) mass->resize(nb_nodes, 0.); if (velocity) velocity->resize(nb_nodes, 0.); if (acceleration) acceleration->resize(nb_nodes, 0.); if (external_force) external_force->resize(nb_nodes, 0.); if (internal_force) internal_force->resize(nb_nodes, 0.); if (blocked_dofs) blocked_dofs->resize(nb_nodes, 0.); if (current_position) current_position->resize(nb_nodes, 0.); if (previous_displacement) previous_displacement->resize(nb_nodes, 0.); if (displacement_increment) displacement_increment->resize(nb_nodes, 0.); for (auto & material : materials) { material->onNodesAdded(nodes_list, event); } if(this->getDOFManager().hasMatrix("M")) { this->assembleMass(nodes_list); } if(this->getDOFManager().hasLumpedMatrix("M")) { this->assembleMassLumped(nodes_list); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesRemoved(__attribute__((unused)) const Array<UInt> & element_list, const Array<UInt> & new_numbering, __attribute__((unused)) const RemovedNodesEvent & event) { if (displacement) { mesh.removeNodesFromArray(*displacement, new_numbering); ++displacement_release; } if (mass) mesh.removeNodesFromArray(*mass, new_numbering); if (velocity) mesh.removeNodesFromArray(*velocity, new_numbering); if (acceleration) mesh.removeNodesFromArray(*acceleration, new_numbering); if (internal_force) mesh.removeNodesFromArray(*internal_force, new_numbering); if (external_force) mesh.removeNodesFromArray(*external_force, new_numbering); if (blocked_dofs) mesh.removeNodesFromArray(*blocked_dofs, new_numbering); // if (increment_acceleration) // mesh.removeNodesFromArray(*increment_acceleration, new_numbering); if (displacement_increment) mesh.removeNodesFromArray(*displacement_increment, new_numbering); if (previous_displacement) mesh.removeNodesFromArray(*previous_displacement, new_numbering); // if (method != _explicit_lumped_mass) { // this->initSolver(); // } } /* -------------------------------------------------------------------------- */ 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 : " << Model::spatial_dimension << std::endl; stream << space << " + fem [" << std::endl; getFEEngine().printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + nodals information [" << std::endl; displacement->printself(stream, indent + 2); mass->printself(stream, indent + 2); velocity->printself(stream, indent + 2); acceleration->printself(stream, indent + 2); external_force->printself(stream, indent + 2); internal_force->printself(stream, indent + 2); blocked_dofs->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + material information [" << std::endl; material_index.printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + materials [" << std::endl; for (auto & material : materials) { material->printself(stream, indent + 1); } stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::insertIntegrationPointsInNeighborhoods( const GhostType & ghost_type) { for (auto & mat : materials) { try { ElementTypeMapArray<Real> quadrature_points_coordinates( "quadrature_points_coordinates_tmp_nl", this->id, this->memory_id); quadrature_points_coordinates.initialize( this->getFEEngine(), _spatial_dimension = Model::spatial_dimension, _nb_component = Model::spatial_dimension, _ghost_type = ghost_type); for (auto & type : quadrature_points_coordinates.elementTypes( Model::spatial_dimension, ghost_type)) { this->getFEEngine().computeIntegrationPointsCoordinates( quadrature_points_coordinates(type, ghost_type), type, ghost_type); } auto & mat_non_local = dynamic_cast<MaterialNonLocalInterface &>(*mat); mat_non_local.insertIntegrationPointsInNeighborhoods( ghost_type, quadrature_points_coordinates); } catch (std::bad_cast &) { } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeNonLocalStresses( const GhostType & ghost_type) { for (auto & mat : materials) { try { auto & mat_non_local = dynamic_cast<MaterialNonLocalInterface &>(*mat); mat_non_local.computeNonLocalStresses(ghost_type); } catch (std::bad_cast &) { } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateLocalInternal( ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind) { const ID field_name = internal_flat.getName(); for (auto & material : materials) { if (material->isInternal<Real>(field_name, kind)) material->flattenInternal(field_name, internal_flat, ghost_type, kind); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateNonLocalInternal( ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind) { const ID field_name = internal_flat.getName(); for (auto & mat : materials) { try { auto & mat_non_local = dynamic_cast<MaterialNonLocalInterface &>(*mat); mat_non_local.updateNonLocalInternals(internal_flat, field_name, ghost_type, kind); } catch (std::bad_cast &) { } } } /* -------------------------------------------------------------------------- */ +FEEngine & SolidMechanicsModel::getFEEngineBoundary(const ID & name) { + return dynamic_cast<FEEngine &>( + getFEEngineClassBoundary<MyFEEngineType>(name)); +} } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh index 4da53ad32..03ce18746 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,821 +1,826 @@ /** * @file solid_mechanics_model.hh * * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> * @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch> * @author Nicolas Richart <nicolas.richart@epfl.ch> * * @date creation: Tue Jul 27 2010 * @date last modification: Tue Jan 19 2016 * * @brief Model of Solid Mechanics * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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 "boundary_condition.hh" #include "data_accessor.hh" #include "model.hh" #include "non_local_manager.hh" #include "solid_mechanics_model_event_handler.hh" +#include "fe_engine.hh" /* -------------------------------------------------------------------------- */ -#include "integrator_gauss.hh" -#include "shape_lagrange.hh" +//#include "integrator_gauss.hh" +//#include "shape_lagrange.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_MECHANICS_MODEL_HH__ #define __AKANTU_SOLID_MECHANICS_MODEL_HH__ namespace akantu { class Material; class MaterialSelector; class DumperIOHelper; class NonLocalManager; +template <ElementKind kind, class IntegrationOrderFunctor> +class IntegratorGauss; +template <ElementKind kind> class ShapeLagrange; } // namespace akantu /* -------------------------------------------------------------------------- */ namespace akantu { struct SolidMechanicsModelOptions : public ModelOptions { explicit SolidMechanicsModelOptions( AnalysisMethod analysis_method = _explicit_lumped_mass); template <typename... pack> SolidMechanicsModelOptions(use_named_args_t, pack &&... _pack); AnalysisMethod analysis_method; - //bool no_init_materials; + // bool no_init_materials; }; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ class SolidMechanicsModel : public Model, public DataAccessor<Element>, public DataAccessor<UInt>, public BoundaryCondition<SolidMechanicsModel>, public NonLocalManagerCallback, public MeshEventHandler, public EventHandlerManager<SolidMechanicsModelEventHandler> { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: class NewMaterialElementsEvent : public NewElementsEvent { public: AKANTU_GET_MACRO_NOT_CONST(MaterialList, material, Array<UInt> &); AKANTU_GET_MACRO(MaterialList, material, const Array<UInt> &); protected: Array<UInt> material; }; using MyFEEngineType = FEEngineTemplate<IntegratorGauss, ShapeLagrange>; protected: using EventManager = EventHandlerManager<SolidMechanicsModelEventHandler>; public: SolidMechanicsModel(Mesh & mesh, UInt spatial_dimension = _all_dimensions, const ID & id = "solid_mechanics_model", const MemoryID & memory_id = 0); virtual ~SolidMechanicsModel(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: template <typename P, typename T, typename... pack> void initFull(named_argument::param_t<P, T &&> && first, pack &&... _pack) { this->initFull(SolidMechanicsModelOptions{use_named_args, first, _pack...}); } /// initialize completely the model void initFull( const ModelOptions & options = SolidMechanicsModelOptions()) override; /// initialize the fem object needed for boundary conditions void initFEEngineBoundary(); /// register the tags associated with the parallel synchronizer // virtual void initParallel(MeshPartition * partition, // DataAccessor<Element> * data_accessor = NULL); /// allocate all vectors virtual void initArrays(); /// allocate all vectors // void initArraysPreviousDisplacment(); /// initialize all internal arrays for materials virtual void initMaterials(); /// initialize the model void initModel() override; /// init PBC synchronizer // void initPBC(); /// initialize a new solver and sets it as the default one to use void initNewSolver(const AnalysisMethod & method); /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; protected: /// allocate an array if needed template <typename T> void allocNodalField(Array<T> *& array, const ID & name); /* ------------------------------------------------------------------------ */ /* PBC */ /* ------------------------------------------------------------------------ */ public: /// change the equation number for proper assembly when using PBC // void changeEquationNumberforPBC(std::map <UInt, UInt> & pbc_pair); /// synchronize Residual for output // void synchronizeResidual(); protected: /// register PBC synchronizer // void registerPBCSynchronizer(); /* ------------------------------------------------------------------------ */ /* Solver interface */ /* ------------------------------------------------------------------------ */ public: /// assembles the stiffness matrix, virtual void assembleStiffnessMatrix(); /// assembles the internal forces in the array internal_forces virtual void assembleInternalForces(); protected: /// callback for the solver, this adds f_{ext} - f_{int} to the residual void assembleResidual() override; /// callback for the solver, this assembles the stiffness matrix void assembleJacobian() override; /// callback for the solver, this is called at beginning of solve void predictor() override; /// callback for the solver, this is called at end of solve void corrector() override; /// Callback for the model to instantiate the matricees when needed void initSolver(TimeStepSolverType time_step_solver_type, NonLinearSolverType non_linear_solver_type) override; protected: /* ------------------------------------------------------------------------ */ TimeStepSolverType getDefaultSolverType() const override; /* ------------------------------------------------------------------------ */ ModelSolverOptions getDefaultSolverOptions(const TimeStepSolverType & type) const override; /* ------------------------------------------------------------------------ */ /* Explicit */ /* ------------------------------------------------------------------------ */ // public: // /// initialize the stuff for the explicit scheme // void initExplicit(AnalysisMethod analysis_method = // _explicit_lumped_mass); public: bool isDefaultSolverExplicit() { return method == _explicit_lumped_mass || method == _explicit_consistent_mass; } // /// initialize the array needed by updateResidual (residual, // current_position) // void initializeUpdateResidualData(); protected: /// update the current position vector void updateCurrentPosition(); // /// assemble the residual for the explicit scheme // virtual void updateResidual(bool need_initialize = true); // /** // * \brief compute the acceleration from the residual // * this function is the explicit equivalent to solveDynamic in implicit // * In the case of lumped mass just divide the residual by the mass // * In the case of not lumped mass call // solveDynamic<_acceleration_corrector> // */ // void updateAcceleration(); /// Update the increment of displacement // void updateIncrement(); // /// Copy the actuel displacement into previous displacement // void updatePreviousDisplacement(); // /// Save stress and strain through EventManager // void saveStressAndStrainBeforeDamage(); // /// Update energies through EventManager // void updateEnergiesAfterDamage(); // /// Solve the system @f[ A x = \alpha b @f] with A a lumped matrix // void solveLumped(Array<Real> & x, const Array<Real> & A, // const Array<Real> & b, const Array<bool> & blocked_dofs, // Real alpha); // /// explicit integration predictor // void explicitPred(); // /// explicit integration corrector // void explicitCorr(); // public: // void solveStep(); // /* // ------------------------------------------------------------------------ // */ // /* Implicit */ // /* // ------------------------------------------------------------------------ // */ // public: // /// initialize the solver and the jacobian_matrix (called by // initImplicit) // void initSolver(); // /// initialize the stuff for the implicit solver // void initImplicit(bool dynamic = false); // /// solve Ma = f to get the initial acceleration // void initialAcceleration(); // /// assemble the stiffness matrix // void assembleStiffnessMatrix(); // public: // /** // * solve a step (predictor + convergence loop + corrector) using the // * the given convergence method (see akantu::SolveConvergenceMethod) // * and the given convergence criteria (see // * akantu::SolveConvergenceCriteria) // **/ // template <SolveConvergenceMethod method, SolveConvergenceCriteria // criteria> // bool solveStep(Real tolerance, UInt max_iteration = 100); // template <SolveConvergenceMethod method, SolveConvergenceCriteria // criteria> // bool solveStep(Real tolerance, Real & error, UInt max_iteration = 100, // bool do_not_factorize = false); // public: // /** // * solve Ku = f using the the given convergence method (see // * akantu::SolveConvergenceMethod) and the given convergence // * criteria (see akantu::SolveConvergenceCriteria) // **/ // template <SolveConvergenceMethod cmethod, SolveConvergenceCriteria // criteria> // bool solveStatic(Real tolerance, UInt max_iteration, // bool do_not_factorize = false); // /// create and return the velocity damping matrix // SparseMatrix & initVelocityDampingMatrix(); // /// implicit time integration predictor // void implicitPred(); // /// implicit time integration corrector // void implicitCorr(); // /// compute the Cauchy stress on user demand. // void computeCauchyStresses(); // // /// compute A and solve @f[ A\delta u = f_ext - f_int @f] // // template <NewmarkBeta::IntegrationSchemeCorrectorType type> // // void solve(Array<Real> &increment, Real block_val = 1., // // bool need_factorize = true, bool has_profile_changed = // false); // protected: // /// finish the computation of residual to solve in increment // void updateResidualInternal(); // /// compute the support reaction and store it in force // void updateSupportReaction(); // private: // /// re-initialize the J matrix (to use if the profile of K changed) // void initJacobianMatrix(); /* ------------------------------------------------------------------------ */ // public: /// Update the stresses for the computation of the residual of the Stiffness /// matrix in the case of finite deformation // void computeStresses(); /// synchronize the ghost element boundaries values // void synchronizeBoundaries(); /* ------------------------------------------------------------------------ */ /* Materials (solid_mechanics_model_material.cc) */ /* ------------------------------------------------------------------------ */ public: /// registers all the custom materials of a given type present in the input /// file // template <typename M> void registerNewCustomMaterials(const ID & mat_type); /// register an empty material of a given type Material & registerNewMaterial(const ID & mat_name, const ID & mat_type, const ID & opt_param); // /// Use a UIntData in the mesh to specify the material to use per // element void setMaterialIDsFromIntData(const std::string & data_name); /// reassigns materials depending on the material selector virtual void reassignMaterial(); /// apply a constant eigen_grad_u on all quadrature points of a given material virtual void applyEigenGradU(const Matrix<Real> & prescribed_eigen_grad_u, const ID & material_name, const GhostType ghost_type = _not_ghost); protected: /// register a material in the dynamic database Material & registerNewMaterial(const ParserSection & mat_section); /// read the material files to instantiate all the materials void instantiateMaterials(); /// set the element_id_by_material and add the elements to the good materials void assignMaterialToElements(const ElementTypeMapArray<UInt> * filter = NULL); /// reinitialize dof_synchronizer and solver (either in implicit or /// explicit) when cohesive elements are inserted // void reinitializeSolver(); /* ------------------------------------------------------------------------ */ /* Mass (solid_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the mass matrix for consistent mass resolutions void assembleMass(); protected: /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumped(GhostType ghost_type); /// assemble the mass matrix for either _ghost or _not_ghost elements void assembleMass(GhostType ghost_type); /// assemble the lumped mass only for a given list of nodes void assembleMassLumped(const Array<UInt> & node_list); /// assemble the lumped mass only for a given list of nodes void assembleMass(const Array<UInt> & node_list); /// fill a vector of rho void computeRho(Array<Real> & rho, ElementType type, GhostType ghost_type); /// compute the kinetic energy Real getKineticEnergy(); Real getKineticEnergy(const ElementType & type, UInt index); /// compute the external work (for impose displacement, the velocity should be /// given too) Real getExternalWork(); /* ------------------------------------------------------------------------ */ /* NonLocalManager inherited members */ /* ------------------------------------------------------------------------ */ protected: void updateDataForNonLocalCriterion(ElementTypeMapReal & criterion) override; void computeNonLocalStresses(const GhostType & ghost_type) override; void insertIntegrationPointsInNeighborhoods(const GhostType & ghost_type) override; /// update the values of the non local internal void updateLocalInternal(ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind) override; /// copy the results of the averaging in the materials void updateNonLocalInternal(ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind) override; /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: inline UInt getNbData(const Array<Element> & elements, const SynchronizationTag & tag) const override; inline void packData(CommunicationBuffer & buffer, const Array<Element> & elements, const SynchronizationTag & tag) const override; inline void unpackData(CommunicationBuffer & buffer, const Array<Element> & elements, const SynchronizationTag & tag) override; inline UInt getNbData(const Array<UInt> & dofs, const SynchronizationTag & tag) const override; inline void packData(CommunicationBuffer & buffer, const Array<UInt> & dofs, const SynchronizationTag & tag) const override; inline void unpackData(CommunicationBuffer & buffer, const Array<UInt> & dofs, const SynchronizationTag & tag) override; protected: inline void splitElementByMaterial(const Array<Element> & elements, Array<Element> * elements_per_mat) const; /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: void onNodesAdded(const Array<UInt> & nodes_list, const NewNodesEvent & event) override; void onNodesRemoved(const Array<UInt> & element_list, const Array<UInt> & new_numbering, const RemovedNodesEvent & event) override; void onElementsAdded(const Array<Element> & nodes_list, const NewElementsEvent & event) override; void onElementsRemoved(const Array<Element> & element_list, const ElementTypeMapArray<UInt> & new_numbering, const RemovedElementsEvent & event) override; void onElementsChanged(const Array<Element> &, const Array<Element> &, const ElementTypeMapArray<UInt> &, const ChangedElementsEvent &) override{}; /* ------------------------------------------------------------------------ */ /* Dumpable interface (kept for convenience) and dumper relative functions */ /* ------------------------------------------------------------------------ */ public: virtual void onDump(); //! decide wether a field is a material internal or not bool isInternal(const std::string & field_name, const ElementKind & element_kind); #ifndef SWIG //! give the amount of data per element virtual ElementTypeMap<UInt> getInternalDataPerElem(const std::string & field_name, const ElementKind & kind); //! flatten a given material internal field ElementTypeMapArray<Real> & flattenInternal(const std::string & field_name, const ElementKind & kind, const GhostType ghost_type = _not_ghost); //! flatten all the registered material internals void flattenAllRegisteredInternals(const ElementKind & kind); #endif dumper::Field * createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) override; dumper::Field * createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) override; dumper::Field * createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind) override; virtual void dump(const std::string & dumper_name); virtual void dump(const std::string & dumper_name, UInt step); virtual void dump(const std::string & dumper_name, Real time, UInt step); void dump() override; virtual void dump(UInt step); virtual void dump(Real time, UInt step); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, Model::spatial_dimension, UInt); /// get the current value of the time step // AKANTU_GET_MACRO(TimeStep, time_step, Real); /// set the value of the time step void setTimeStep(Real time_step, const ID & solver_id = "") override; /// void setTimeStep(Real time_step); /// return the of iterations done in the last solveStep // AKANTU_GET_MACRO(NumberIter, n_iter, UInt); /// get the value of the conversion from forces/ mass to acceleration AKANTU_GET_MACRO(F_M2A, f_m2a, Real); /// set the value of the conversion from forces/ mass to acceleration AKANTU_SET_MACRO(F_M2A, f_m2a, Real); /// get the SolidMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement, Array<Real> &); /// get the SolidMechanicsModel::previous_displacement vector AKANTU_GET_MACRO(PreviousDisplacement, *previous_displacement, Array<Real> &); /// get the SolidMechanicsModel::current_position vector \warn only consistent /// after a call to SolidMechanicsModel::updateCurrentPosition const Array<Real> & getCurrentPosition(); /// get the SolidMechanicsModel::increment vector \warn only consistent if /// SolidMechanicsModel::setIncrementFlagOn has been called before AKANTU_GET_MACRO(Increment, *displacement_increment, Array<Real> &); /// get the lumped SolidMechanicsModel::mass vector AKANTU_GET_MACRO(Mass, *mass, Array<Real> &); /// get the SolidMechanicsModel::velocity vector AKANTU_GET_MACRO(Velocity, *velocity, Array<Real> &); /// get the SolidMechanicsModel::acceleration vector, updated by /// SolidMechanicsModel::updateAcceleration AKANTU_GET_MACRO(Acceleration, *acceleration, Array<Real> &); /// get the SolidMechanicsModel::force vector (external forces) AKANTU_GET_MACRO(Force, *external_force, Array<Real> &); /// get the SolidMechanicsModel::internal_force vector (internal forces) AKANTU_GET_MACRO(InternalForce, *internal_force, Array<Real> &); /// get the SolidMechanicsModel::blocked_dofs vector AKANTU_GET_MACRO(BlockedDOFs, *blocked_dofs, Array<bool> &); /// get the SolidMechnicsModel::incrementAcceleration vector // AKANTU_GET_MACRO(IncrementAcceleration, *increment_acceleration, // Array<Real> &); /// get the value of the SolidMechanicsModel::increment_flag AKANTU_GET_MACRO(IncrementFlag, increment_flag, bool); /// get a particular material (by material index) inline Material & getMaterial(UInt mat_index); /// get a particular material (by material index) inline const Material & getMaterial(UInt mat_index) const; /// get a particular material (by material name) inline Material & getMaterial(const std::string & name); /// get a particular material (by material name) inline const Material & getMaterial(const std::string & name) const; /// get a particular material id from is name inline UInt getMaterialIndex(const std::string & name) const; /// give the number of materials inline UInt getNbMaterials() const { return materials.size(); } inline void setMaterialSelector(MaterialSelector & selector); /// give the material internal index from its id Int getInternalIndexFromID(const ID & id) const; /// compute the stable time step Real getStableTimeStep(); /// get the energies Real getEnergy(const std::string & energy_id); /// compute the energy for energy Real getEnergy(const std::string & energy_id, const ElementType & type, UInt index); /** * @brief set the SolidMechanicsModel::increment_flag to on, the activate the * update of the SolidMechanicsModel::increment vector */ // void setIncrementFlagOn(); // /// get the stiffness matrix // AKANTU_GET_MACRO(StiffnessMatrix, *stiffness_matrix, SparseMatrix &); // /// get the global jacobian matrix of the system // AKANTU_GET_MACRO(GlobalJacobianMatrix, *jacobian_matrix, const SparseMatrix // &); // /// get the mass matrix // AKANTU_GET_MACRO(MassMatrix, *mass_matrix, SparseMatrix &); // /// get the velocity damping matrix // AKANTU_GET_MACRO(VelocityDampingMatrix, *velocity_damping_matrix, // SparseMatrix &); /// get the FEEngine object to integrate or interpolate on the boundary - inline FEEngine & getFEEngineBoundary(const ID & name = ""); + FEEngine & getFEEngineBoundary(const ID & name = ""); // /// get integrator // AKANTU_GET_MACRO(Integrator, *integrator, const IntegrationScheme2ndOrder // &); // /// get synchronizer // AKANTU_GET_MACRO(Synchronizer, *synch_parallel, const ElementSynchronizer // &); AKANTU_GET_MACRO(MaterialByElement, material_index, const ElementTypeMapArray<UInt> &); /// vectors containing local material element index for each global element /// index AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialByElement, material_index, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialByElement, material_index, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(MaterialLocalNumbering, material_local_numbering, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(MaterialLocalNumbering, material_local_numbering, UInt); /// Get the type of analysis method used AKANTU_GET_MACRO(AnalysisMethod, method, AnalysisMethod); /// Access the non_local_manager interface AKANTU_GET_MACRO(NonLocalManager, *non_local_manager, NonLocalManager &); + protected: friend class Material; protected: /// compute the stable time step Real getStableTimeStep(const GhostType & ghost_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// number of iterations // UInt n_iter; /// time step // Real time_step; /// conversion coefficient form force/mass to acceleration Real f_m2a; /// displacements array Array<Real> * displacement; UInt displacement_release{0}; /// displacements array at the previous time step (used in finite deformation) Array<Real> * previous_displacement; /// increment of displacement Array<Real> * displacement_increment; /// lumped mass array Array<Real> * mass; /// velocities array Array<Real> * velocity; /// accelerations array Array<Real> * acceleration; /// accelerations array // Array<Real> * increment_acceleration; /// external forces array Array<Real> * external_force; /// internal forces array Array<Real> * internal_force; /// array specifing if a degree of freedom is blocked or not Array<bool> * blocked_dofs; /// array of current position used during update residual Array<Real> * current_position; UInt current_position_release{0}; /// mass matrix SparseMatrix * mass_matrix; /// velocity damping matrix SparseMatrix * velocity_damping_matrix; /// stiffness matrix SparseMatrix * stiffness_matrix; /// jacobian matrix @f[A = cM + dD + K@f] with @f[c = \frac{1}{\beta \Delta /// t^2}, d = \frac{\gamma}{\beta \Delta t} @f] SparseMatrix * jacobian_matrix; /// Arrays containing the material index for each element ElementTypeMapArray<UInt> material_index; /// Arrays containing the position in the element filter of the material /// (material's local numbering) ElementTypeMapArray<UInt> material_local_numbering; #ifdef SWIGPYTHON public: #endif /// list of used materials std::vector<std::unique_ptr<Material>> materials; /// mapping between material name and material internal id std::map<std::string, UInt> materials_names_to_id; #ifdef SWIGPYTHON protected: #endif /// class defining of to choose a material MaterialSelector * material_selector; /// define if it is the default selector or not bool is_default_material_selector; // /// integration scheme of second order used // IntegrationScheme2ndOrder *integrator; /// flag defining if the increment must be computed or not bool increment_flag; /// analysis method check the list in akantu::AnalysisMethod AnalysisMethod method; /// tells if the material are instantiated bool are_materials_instantiated; typedef std::map<std::pair<std::string, ElementKind>, ElementTypeMapArray<Real> *> flatten_internal_map; /// map a registered internals to be flattened for dump purposes flatten_internal_map registered_internals; /// non local manager std::unique_ptr<NonLocalManager> non_local_manager; /// pointer to the pbc synchronizer // PBCSynchronizer * pbc_synch; }; /* -------------------------------------------------------------------------- */ namespace BC { namespace Neumann { typedef FromHigherDim FromStress; typedef FromSameDim FromTraction; } // namespace Neumann } // namespace BC /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const SolidMechanicsModel & _this) { _this.printself(stream); return stream; } } // namespace akantu /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "parser.hh" #include "solid_mechanics_model_inline_impl.cc" #include "solid_mechanics_model_tmpl.hh" #include "material_selector_tmpl.hh" #endif /* __AKANTU_SOLID_MECHANICS_MODEL_HH__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.cc index adbb273a4..979a6c0f4 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.cc @@ -1,615 +1,615 @@ /** * @file fragment_manager.cc * * @author Marco Vocialta <marco.vocialta@epfl.ch> * * @date creation: Thu Jan 23 2014 * @date last modification: Mon Dec 14 2015 * * @brief Group manager to handle fragments * * @section LICENSE * * Copyright (©) 2014, 2015 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 "aka_iterators.hh" #include "material_cohesive.hh" +#include "mesh_iterators.hh" #include "solid_mechanics_model_cohesive.hh" #include "static_communicator.hh" /* -------------------------------------------------------------------------- */ #include <algorithm> #include <functional> #include <numeric> /* -------------------------------------------------------------------------- */ namespace 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"), principal_directions( 0, model.getSpatialDimension() * model.getSpatialDimension(), "principal_directions"), quad_coordinates("quad_coordinates", id), mass_density("mass_density", id), nb_elements_per_fragment(0, 1, "nb_elements_per_fragment"), dump_data(dump_data) { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); /// compute quadrature points' coordinates quad_coordinates.initialize(mesh, _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension, _ghost_type = _not_ghost); // mesh.initElementTypeMapArray(quad_coordinates, spatial_dimension, // spatial_dimension, _not_ghost); model.getFEEngine().interpolateOnIntegrationPoints(model.getMesh().getNodes(), quad_coordinates); /// store mass density per quadrature point mass_density.initialize(mesh, _spatial_dimension = spatial_dimension, _ghost_type = _not_ghost); // mesh.initElementTypeMapArray(mass_density, 1, spatial_dimension, // _not_ghost); storeMassDensityPerIntegrationPoint(); 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> & mat_indexes = model.getMaterialByElement(el.type, el.ghost_type); const Array<UInt> & mat_loc_num = model.getMaterialLocalNumbering(el.type, el.ghost_type); const MaterialCohesive & mat = static_cast<const MaterialCohesive &>( model.getMaterial(mat_indexes(el.element))); UInt el_index = mat_loc_num(el.element); UInt nb_quad_per_element = model.getFEEngine("CohesiveFEEngine") .getNbIntegrationPoints(el.type, el.ghost_type); const Array<Real> & damage_array = mat.getDamage(el.type, el.ghost_type); AKANTU_DEBUG_ASSERT(nb_quad_per_element * el_index < damage_array.getSize(), "This quadrature point is out of range"); const Real * element_damage = damage_array.storage() + nb_quad_per_element * el_index; UInt unbroken_quads = std::count_if( element_damage, element_damage + 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(Real damage_limit) { AKANTU_DEBUG_IN(); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) ElementSynchronizer * cohesive_synchronizer = const_cast<ElementSynchronizer *>(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 auto & mesh_facets = const_cast<Mesh &>(mesh.getMeshFacets()); UInt spatial_dimension = model.getSpatialDimension(); std::string fragment_prefix("fragment"); /// generate fragments global_nb_fragment = - createClusters(spatial_dimension, mesh_facets, fragment_prefix, + createClusters(spatial_dimension, mesh_facets, fragment_prefix, CohesiveElementFilter(model, damage_limit)); 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", true); 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 ElementTypeMapArray<Real> unit_field("unit_field", id); unit_field.initialize(model.getFEEngine(), _nb_component = spatial_dimension, _spatial_dimension = spatial_dimension, _ghost_type = _not_ghost, _default_value = 1.); // mesh.initElementTypeMapArray(unit_field, spatial_dimension, // spatial_dimension, // _not_ghost); // ElementTypeMapArray<Real>::type_iterator it = // unit_field.firstType(spatial_dimension, _not_ghost, _ek_regular); // ElementTypeMapArray<Real>::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.getFEEngine().getNbIntegrationPoints(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 ElementTypeMapArray<Real> velocity_field("velocity_field", id); - velocity_field.initialize(model.getFEEngine(), _nb_component = spatial_dimension, - _spatial_dimension = spatial_dimension, - _ghost_type = _not_ghost); + velocity_field.initialize( + model.getFEEngine(), _nb_component = spatial_dimension, + _spatial_dimension = spatial_dimension, _ghost_type = _not_ghost); // mesh.initElementTypeMapArray(velocity_field, spatial_dimension, // spatial_dimension, _not_ghost); model.getFEEngine().interpolateOnIntegrationPoints(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 ElementTypeMapArray<Real> moments_coords("moments_coords", id); - moments_coords.initialize(model.getFEEngine(), _nb_component = spatial_dimension * spatial_dimension, - _spatial_dimension = spatial_dimension, - _ghost_type = _not_ghost, _default_value = 1.); + moments_coords.initialize(model.getFEEngine(), + _nb_component = + spatial_dimension * spatial_dimension, + _spatial_dimension = spatial_dimension, + _ghost_type = _not_ghost, _default_value = 1.); // mesh.initElementTypeMapArray(moments_coords, // spatial_dimension * spatial_dimension, // spatial_dimension, _not_ghost); // /// resize the by element type // ElementTypeMapArray<Real>::type_iterator it = // moments_coords.firstType(spatial_dimension, _not_ghost, _ek_regular); // ElementTypeMapArray<Real>::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.getFEEngine().getNbIntegrationPoints(type); + // UInt nb_quad_per_element = + // model.getFEEngine().getNbIntegrationPoints(type); // field_array.resize(nb_element * nb_quad_per_element); // } /// compute coordinates Array<Real>::const_vector_iterator 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 ElementTypeMapArray<UInt> & el_list = it->second->getElements(); /// loop over elements of the fragment - for (auto type : el_list.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { + for (auto type : + el_list.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { UInt nb_quad_per_element = model.getFEEngine().getNbIntegrationPoints(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>::matrix_iterator moments_begin = moments_coords_array.begin(spatial_dimension, spatial_dimension); Array<Real>::const_vector_iterator 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); principal_directions.resize(global_nb_fragment); Array<Real>::matrix_iterator integrated_moments_it = integrated_moments.begin(spatial_dimension, spatial_dimension); Array<Real>::vector_iterator inertia_moments_it = inertia_moments.begin(spatial_dimension); Array<Real>::matrix_iterator principal_directions_it = principal_directions.begin(spatial_dimension, spatial_dimension); for (UInt frag = 0; frag < global_nb_fragment; ++frag, ++integrated_moments_it, ++inertia_moments_it, ++principal_directions_it) { integrated_moments_it->eig(*inertia_moments_it, *principal_directions_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::storeMassDensityPerIntegrationPoint() { 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.getFEEngine().getNbIntegrationPoints(type); mass_density_array.resize(nb_element * nb_quad_per_element); const Array<UInt> & mat_indexes = model.getMaterialByElement(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(mat_indexes(el)); for (UInt q = 0; q < nb_quad_per_element; ++q, ++mass_density_it) *mass_density_it = mat.getRho(); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::integrateFieldOnFragments( ElementTypeMapArray<Real> & 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>::vector_iterator 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 ElementTypeMapArray<UInt> & el_list = it->second->getElements(); /// loop over elements of the fragment for (auto type : el_list.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { UInt nb_quad_per_element = model.getFEEngine().getNbIntegrationPoints(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>::matrix_iterator int_array_it = integration_array.begin_reinterpret(nb_quad_per_element, nb_component, nb_element); Array<Real>::matrix_iterator int_array_end = integration_array.end_reinterpret(nb_quad_per_element, nb_component, nb_element); Array<Real>::matrix_iterator field_array_begin = field_array.begin_reinterpret(nb_quad_per_element, nb_component, field_array.getSize() / nb_quad_per_element); Array<Real>::const_vector_iterator 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.getFEEngine().integrate(integration_array, integrated_array, nb_component, type, _not_ghost, elements); /// sum over all elements and store the result Vector<Real> output_tmp(output_begin[*fragment_index_it]); output_tmp += std::accumulate(integrated_array.begin(nb_component), integrated_array.end(nb_component), Vector<Real>(nb_component)); } } /// sum output over all processors const StaticCommunicator & comm = mesh.getCommunicator(); comm.allReduce(output, _so_sum); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FragmentManager::computeNbElementsPerFragment() { AKANTU_DEBUG_IN(); UInt spatial_dimension = model.getSpatialDimension(); nb_elements_per_fragment.resize(global_nb_fragment); nb_elements_per_fragment.clear(); 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) { const ElementTypeMapArray<UInt> & el_list = it->second->getElements(); /// loop over elements of the fragment for (auto type : el_list.elementTypes(spatial_dimension, _not_ghost, _ek_regular)) { UInt nb_element = el_list(type).getSize(); nb_elements_per_fragment(*fragment_index_it) += nb_element; } } /// sum values over all processors const auto & comm = mesh.getCommunicator(); comm.allReduce(nb_elements_per_fragment, _so_sum); if (dump_data) createDumpDataArray(nb_elements_per_fragment, "elements per fragment"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template <typename T> void FragmentManager::createDumpDataArray(Array<T> & data, std::string name, bool fragment_index_output) { AKANTU_DEBUG_IN(); if (data.getSize() == 0) return; - 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(); + auto & mesh_not_const = const_cast<Mesh &>(mesh); - auto data_begin = data.begin(nb_component); + auto && spatial_dimension = mesh.getSpatialDimension(); + auto && nb_component = data.getNbComponent(); + auto && data_begin = data.begin(nb_component); + auto fragment_index_it = fragment_index.begin(); /// loop over fragments - for (auto it = element_group_begin(); - it != element_group_end(); ++it, ++fragment_index_it) { - - const auto & fragment = *(it->second); + for (const auto & fragment : ElementGroupsIterable(*this)) { + const auto & fragment_idx = *fragment_index_it; /// loop over cluster types for (auto & type : fragment.elementTypes(spatial_dimension)) { /// init mesh data auto & mesh_data = mesh_not_const.getDataPointer<T>( name, type, _not_ghost, nb_component); auto mesh_data_begin = mesh_data.begin(nb_component); /// fill mesh data - if (fragment_index_output) { - for (const auto & elem : fragment.getElements(type)) { - Vector<T> md_tmp(mesh_data_begin[elem]); - md_tmp(0) = *fragment_index_it; - } - } else { - for (const auto & elem : fragment.getElements(type)) { - Vector<T> md_tmp(mesh_data_begin[elem]); - md_tmp = data_begin[*fragment_index_it]; + for (const auto & elem : fragment.getElements(type)) { + Vector<T> md_tmp = mesh_data_begin[elem]; + if (fragment_index_output) { + md_tmp(0) = fragment_idx; + } else { + md_tmp = data_begin[fragment_idx]; } } } } AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc index efa400cbe..f4fc4eb29 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc @@ -1,788 +1,788 @@ /** * @file solid_mechanics_model_cohesive.cc * * @author Mauro Corrado <mauro.corrado@epfl.ch> * @author Nicolas Richart <nicolas.richart@epfl.ch> * @author Marco Vocialta <marco.vocialta@epfl.ch> * * @date creation: Tue May 08 2012 * @date last modification: Wed Jan 13 2016 * * @brief Solid mechanics model for cohesive elements * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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 "solid_mechanics_model_cohesive.hh" #include "aka_iterators.hh" #include "dumpable_inline_impl.hh" #include "material_cohesive.hh" #include "shape_cohesive.hh" #ifdef AKANTU_USE_IOHELPER #include "dumper_paraview.hh" #endif /* -------------------------------------------------------------------------- */ #include <algorithm> /* -------------------------------------------------------------------------- */ namespace akantu { const SolidMechanicsModelCohesiveOptions default_solid_mechanics_model_cohesive_options(_explicit_lumped_mass, false); /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::SolidMechanicsModelCohesive( Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : SolidMechanicsModel(mesh, dim, id, memory_id), tangents("tangents", id), facet_stress("facet_stress", id), facet_material("facet_material", id) { AKANTU_DEBUG_IN(); inserter = NULL; #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) facet_synchronizer = NULL; facet_stress_synchronizer = NULL; cohesive_element_synchronizer = NULL; global_connectivity = NULL; #endif delete material_selector; material_selector = new DefaultMaterialCohesiveSelector(*this); this->registerEventHandler(*this, _ehp_solid_mechanics_model_cohesive); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper<DumperParaview>("cohesive elements", id); this->mesh.addDumpMeshToDumper("cohesive elements", mesh, Model::spatial_dimension, _not_ghost, _ek_cohesive); #endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModelCohesive::~SolidMechanicsModelCohesive() { AKANTU_DEBUG_IN(); delete inserter; #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) delete cohesive_element_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) this->mesh.getDumper("cohesive elements").setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initFull(const ModelOptions & options) { AKANTU_DEBUG_IN(); const SolidMechanicsModelCohesiveOptions & smmc_options = dynamic_cast<const SolidMechanicsModelCohesiveOptions &>(options); this->is_extrinsic = smmc_options.extrinsic; if (!inserter) inserter = new CohesiveElementInserter(mesh, is_extrinsic, id + ":cohesive_element_inserter"); SolidMechanicsModel::initFull(options); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initMaterials() { AKANTU_DEBUG_IN(); // 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].get()) == nullptr) && 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 if (is_extrinsic) { const Mesh & mesh_facets = inserter->getMeshFacets(); facet_material.initialize(mesh_facets, _spatial_dimension = Model::spatial_dimension - 1); // mesh_facets.initElementTypeMapArray(facet_material, 1, // spatial_dimension - 1); Element element; for (auto ghost_type : ghost_types) { element.ghost_type = ghost_type; for (auto & type : mesh_facets.elementTypes(Model::spatial_dimension - 1, ghost_type)) { element.type = type; Array<UInt> & f_material = facet_material(type, ghost_type); UInt nb_element = mesh_facets.getNbElement(type, ghost_type); 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); } } } SolidMechanicsModel::initMaterials(); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (facet_synchronizer != NULL) inserter->initParallel(facet_synchronizer, cohesive_element_synchronizer); // inserter->initParallel(facet_synchronizer, synch_parallel); #endif initAutomaticInsertion(); } else { // TODO think of something a bit mor consistant than just coding the first // thing that comes in Fabian's head.... typedef ParserSection::const_section_iterator const_section_iterator; std::pair<const_section_iterator, const_section_iterator> sub_sections = this->parser->getSubSections(_st_mesh); if (sub_sections.first != sub_sections.second) { std::string cohesive_surfaces = sub_sections.first->getParameter("cohesive_surfaces"); this->initIntrinsicCohesiveMaterials(cohesive_surfaces); } else { this->initIntrinsicCohesiveMaterials(cohesive_index); } } AKANTU_DEBUG_OUT(); } // namespace akantu /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initIntrinsicCohesiveMaterials( std::string cohesive_surfaces) { AKANTU_DEBUG_IN(); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (facet_synchronizer != NULL) inserter->initParallel(facet_synchronizer, cohesive_element_synchronizer); // inserter->initParallel(facet_synchronizer, synch_parallel); #endif std::istringstream split(cohesive_surfaces); std::string physname; while (std::getline(split, physname, ',')) { AKANTU_DEBUG_INFO( "Pre-inserting cohesive elements along facets from physical surface: " << physname); insertElementsFromMeshData(physname); } synchronizeInsertionData(); SolidMechanicsModel::initMaterials(); if (is_default_material_selector) delete material_selector; material_selector = new MeshDataMaterialCohesiveSelector(*this); // UInt nb_new_elements = inserter->insertElements(); // if (nb_new_elements > 0) { // this->reinitializeSolver(); // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::synchronizeInsertionData() { #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (facet_synchronizer != NULL) { facet_synchronizer->asynchronousSynchronize(*inserter, _gst_ce_groups); facet_synchronizer->waitEndSynchronize(*inserter, _gst_ce_groups); } #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initIntrinsicCohesiveMaterials( UInt cohesive_index) { AKANTU_DEBUG_IN(); for (ghost_type_t::iterator gt = ghost_type_t::begin(); gt != ghost_type_t::end(); ++gt) { Mesh::type_iterator first = mesh.firstType(Model::spatial_dimension, *gt, _ek_cohesive); Mesh::type_iterator last = mesh.lastType(Model::spatial_dimension, *gt, _ek_cohesive); for (; first != last; ++first) { Array<UInt> & mat_indexes = this->material_index(*first, *gt); Array<UInt> & mat_loc_num = this->material_local_numbering(*first, *gt); mat_indexes.set(cohesive_index); mat_loc_num.clear(); } } #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (facet_synchronizer != NULL) inserter->initParallel(facet_synchronizer, cohesive_element_synchronizer); // inserter->initParallel(facet_synchronizer, synch_parallel); #endif 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(); registerFEEngineObject<MyFEEngineCohesiveType>("CohesiveFEEngine", mesh, Model::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(Model::spatial_dimension, type_ghost); Mesh::type_iterator last = mesh.lastType(Model::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 = FEEngine::getCohesiveElementType(type_facet); mesh.addConnectivityType(type_cohesive, type_ghost); } } } AKANTU_DEBUG_ASSERT(type != _not_defined, "No elements in the mesh"); getFEEngine("CohesiveFEEngine").initShapeFunctions(_not_ghost); getFEEngine("CohesiveFEEngine").initShapeFunctions(_ghost); registerFEEngineObject<MyFEEngineFacetType>( "FacetsFEEngine", mesh.getMeshFacets(), Model::spatial_dimension - 1); if (is_extrinsic) { getFEEngine("FacetsFEEngine").initShapeFunctions(_not_ghost); getFEEngine("FacetsFEEngine").initShapeFunctions(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::limitInsertion(BC::Axis axis, Real first_limit, Real second_limit) { AKANTU_DEBUG_IN(); inserter->setLimit(axis, first_limit, second_limit); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::insertIntrinsicElements() { AKANTU_DEBUG_IN(); // UInt nb_new_elements = inserter->insertIntrinsicElements(); // if (nb_new_elements > 0) { // this->reinitializeSolver(); // } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::insertElementsFromMeshData( std::string physname) { AKANTU_DEBUG_IN(); UInt material_index = SolidMechanicsModel::getMaterialIndex(physname); inserter->insertIntrinsicElements(physname, material_index); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initAutomaticInsertion() { AKANTU_DEBUG_IN(); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (facet_stress_synchronizer != NULL) { DataAccessor * data_accessor = this; const ElementTypeMapArray<UInt> & rank_to_element = synch_parallel->getPrankToElement(); facet_stress_synchronizer->updateFacetStressSynchronizer( *inserter, rank_to_element, *data_accessor); } #endif facet_stress.initialize(inserter->getMeshFacets(), _nb_component = 2 * Model::spatial_dimension * Model::spatial_dimension, _spatial_dimension = Model::spatial_dimension - 1); // inserter->getMeshFacets().initElementTypeMapArray( // facet_stress, 2 * spatial_dimension * spatial_dimension, // spatial_dimension - 1); resizeFacetStress(); /// compute normals on facets computeNormals(); 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 ElementTypeMapArray<UInt> & rank_to_element = synch_parallel->getPrankToElement(); facet_stress_synchronizer->updateFacetStressSynchronizer( *inserter, rank_to_element, *data_accessor); } #endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::initStressInterpolation() { Mesh & mesh_facets = inserter->getMeshFacets(); /// compute quadrature points coordinates on facets Array<Real> & position = mesh.getNodes(); ElementTypeMapArray<Real> quad_facets("quad_facets", id); quad_facets.initialize(mesh_facets, _nb_component = Model::spatial_dimension, _spatial_dimension = Model::spatial_dimension - 1); // mesh_facets.initElementTypeMapArray(quad_facets, Model::spatial_dimension, // Model::spatial_dimension - 1); getFEEngine("FacetsFEEngine") .interpolateOnIntegrationPoints(position, quad_facets); /// compute elements quadrature point positions and build /// element-facet quadrature points data structure ElementTypeMapArray<Real> elements_quad_facets("elements_quad_facets", id); elements_quad_facets.initialize( mesh, _nb_component = Model::spatial_dimension, _spatial_dimension = Model::spatial_dimension); // mesh.initElementTypeMapArray(elements_quad_facets, // Model::spatial_dimension, // Model::spatial_dimension); for (auto elem_gt : ghost_types) { for (auto & type : mesh.elementTypes(Model::spatial_dimension, elem_gt)) { 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 = getFEEngine("FacetsFEEngine").getNbIntegrationPoints(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 < Model::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::assembleInternalForces() { AKANTU_DEBUG_IN(); // f_int += f_int_cohe for (auto & material : this->materials) { try { MaterialCohesive & mat = dynamic_cast<MaterialCohesive &>(*material); mat.computeTraction(_not_ghost); } catch (std::bad_cast & bce) { } } SolidMechanicsModel::assembleInternalForces(); if (isDefaultSolverExplicit()) { for (auto & material : materials) { try { MaterialCohesive & mat = dynamic_cast<MaterialCohesive &>(*material); mat.computeEnergies(); } catch (std::bad_cast & bce) { } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::computeNormals() { AKANTU_DEBUG_IN(); Mesh & mesh_facets = this->inserter->getMeshFacets(); this->getFEEngine("FacetsFEEngine") .computeNormalsOnIntegrationPoints(_not_ghost); /** * @todo store tangents while computing normals instead of * recomputing them as follows: */ /* ------------------------------------------------------------------------ */ UInt tangent_components = Model::spatial_dimension * (Model::spatial_dimension - 1); tangents.initialize(mesh_facets, _nb_component = tangent_components, _spatial_dimension = Model::spatial_dimension - 1); // mesh_facets.initElementTypeMapArray(tangents, tangent_components, // Model::spatial_dimension - 1); for (auto facet_type : mesh_facets.elementTypes(Model::spatial_dimension - 1)) { const Array<Real> & normals = this->getFEEngine("FacetsFEEngine") .getNormalsOnIntegrationPoints(facet_type); Array<Real> & tangents = this->tangents(facet_type); Math::compute_tangents(normals, tangents); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::interpolateStress() { ElementTypeMapArray<Real> by_elem_result("temporary_stress_by_facets", id); for (auto & material : materials) { try { MaterialCohesive & mat __attribute__((unused)) = dynamic_cast<MaterialCohesive &>(*material); } catch (std::bad_cast &) { /// interpolate stress on facet quadrature points positions material->interpolateStressOnFacets(facet_stress, by_elem_result); } } #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData( debug::_dm_model_cohesive, "Interpolated stresses before", facet_stress); #endif this->synchronize(_gst_smmc_facets_stress); #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData(debug::_dm_model_cohesive, "Interpolated stresses", facet_stress); #endif } /* -------------------------------------------------------------------------- */ UInt SolidMechanicsModelCohesive::checkCohesiveStress() { interpolateStress(); for (auto & mat : materials) { try { MaterialCohesive & mat_cohesive = dynamic_cast<MaterialCohesive &>(*mat); /// check which not ghost cohesive elements are to be created mat_cohesive.checkInsertion(); } catch (std::bad_cast &) { } } /// communicate data among processors this->synchronize(_gst_smmc_facets); /// insert cohesive elements UInt nb_new_elements = inserter->insertElements(); // if (nb_new_elements > 0) { // this->reinitializeSolver(); // } return nb_new_elements; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onElementsAdded( const Array<Element> & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) updateCohesiveSynchronizers(); #endif SolidMechanicsModel::onElementsAdded(element_list, event); #if defined(AKANTU_PARALLEL_COHESIVE_ELEMENT) if (cohesive_element_synchronizer != NULL) cohesive_element_synchronizer->computeAllBufferSizes(*this); #endif if (is_extrinsic) resizeFacetStress(); /// if (method != _explicit_lumped_mass) { /// this->initSolver(); /// } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onNodesAdded(const Array<UInt> & new_nodes, const NewNodesEvent & event) { AKANTU_DEBUG_IN(); // 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(new_nodes, event); UInt new_node, old_node; try { const auto & cohesive_event = dynamic_cast<const CohesiveNewNodesEvent &>(event); const auto & old_nodes = cohesive_event.getOldNodesList(); auto copy = [this, &new_node, &old_node](auto & arr) { for (UInt s = 0; s < spatial_dimension; ++s) { arr(new_node, s) = arr(old_node, s); } }; - for (auto pair : zip(new_nodes, old_nodes)) { + for (auto && pair : zip(new_nodes, old_nodes)) { std::tie(new_node, old_node) = pair; copy(*displacement); copy(*velocity); copy(*acceleration); copy(*blocked_dofs); if (current_position) copy(*current_position); if (previous_displacement) copy(*previous_displacement); } // if (this->getDOFManager().hasMatrix("M")) { // this->assembleMass(old_nodes); // } // if (this->getDOFManager().hasLumpedMatrix("M")) { // this->assembleMassLumped(old_nodes); // } } catch (std::bad_cast &) { } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onEndSolveStep(const AnalysisMethod &) { AKANTU_DEBUG_IN(); /* * This is required because the Cauchy stress is the stress measure that * is used to check the insertion of cohesive elements */ for (auto & mat : materials) { if (mat->isFiniteDeformation()) mat->computeAllCauchyStresses(_not_ghost); } 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(); Mesh & mesh_facets = inserter->getMeshFacets(); 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(Model::spatial_dimension - 1, ghost_type); Mesh::type_iterator end = mesh_facets.lastType(Model::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 = getFEEngine("FacetsFEEngine") .getNbIntegrationPoints(type, ghost_type); UInt new_size = nb_facet * nb_quadrature_points; facet_stress(type, ghost_type).resize(new_size); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::addDumpGroupFieldToDumper( const std::string & dumper_name, const std::string & field_id, const std::string & group_name, const ElementKind & element_kind, bool padding_flag) { AKANTU_DEBUG_IN(); UInt spatial_dimension = Model::spatial_dimension; ElementKind _element_kind = element_kind; if (dumper_name == "cohesive elements") { _element_kind = _ek_cohesive; } else if (dumper_name == "facets") { spatial_dimension = Model::spatial_dimension - 1; } SolidMechanicsModel::addDumpGroupFieldToDumper(dumper_name, field_id, group_name, spatial_dimension, _element_kind, padding_flag); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModelCohesive::onDump() { this->flattenAllRegisteredInternals(_ek_cohesive); SolidMechanicsModel::onDump(); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_inline_impl.cc b/src/model/solid_mechanics/solid_mechanics_model_inline_impl.cc index 4e25fafbd..24f214852 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_inline_impl.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_inline_impl.cc @@ -1,426 +1,421 @@ /** * @file solid_mechanics_model_inline_impl.cc * * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> * @author Daniel Pino Muñoz <daniel.pinomunoz@epfl.ch> * @author Nicolas Richart <nicolas.richart@epfl.ch> * * @date creation: Wed Aug 04 2010 * @date last modification: Wed Nov 18 2015 * * @brief Implementation of the inline functions of the SolidMechanicsModel * class * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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_named_argument.hh" #include "material_selector.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_MECHANICS_MODEL_INLINE_IMPL_CC__ #define __AKANTU_SOLID_MECHANICS_MODEL_INLINE_IMPL_CC__ namespace akantu { /* -------------------------------------------------------------------------- */ inline SolidMechanicsModelOptions::SolidMechanicsModelOptions( AnalysisMethod analysis_method) : analysis_method(analysis_method) {} /* -------------------------------------------------------------------------- */ template <typename... pack> SolidMechanicsModelOptions::SolidMechanicsModelOptions(use_named_args_t, pack &&... _pack) : SolidMechanicsModelOptions( OPTIONAL_NAMED_ARG(analysis_method, _explicit_lumped_mass)) {} /* -------------------------------------------------------------------------- */ inline Material & SolidMechanicsModel::getMaterial(UInt mat_index) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(mat_index < materials.size(), "The model " << id << " has no material no " << mat_index); AKANTU_DEBUG_OUT(); return *materials[mat_index]; } /* -------------------------------------------------------------------------- */ inline const Material & SolidMechanicsModel::getMaterial(UInt mat_index) const { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(mat_index < materials.size(), "The model " << id << " has no material no " << mat_index); AKANTU_DEBUG_OUT(); return *materials[mat_index]; } /* -------------------------------------------------------------------------- */ inline Material & SolidMechanicsModel::getMaterial(const std::string & name) { AKANTU_DEBUG_IN(); std::map<std::string, UInt>::const_iterator it = materials_names_to_id.find(name); AKANTU_DEBUG_ASSERT(it != materials_names_to_id.end(), "The model " << id << " has no material named " << name); AKANTU_DEBUG_OUT(); return *materials[it->second]; } /* -------------------------------------------------------------------------- */ inline UInt SolidMechanicsModel::getMaterialIndex(const std::string & name) const { AKANTU_DEBUG_IN(); std::map<std::string, UInt>::const_iterator it = materials_names_to_id.find(name); AKANTU_DEBUG_ASSERT(it != materials_names_to_id.end(), "The model " << id << " has no material named " << name); AKANTU_DEBUG_OUT(); return it->second; } /* -------------------------------------------------------------------------- */ inline const Material & SolidMechanicsModel::getMaterial(const std::string & name) const { AKANTU_DEBUG_IN(); std::map<std::string, UInt>::const_iterator it = materials_names_to_id.find(name); AKANTU_DEBUG_ASSERT(it != materials_names_to_id.end(), "The model " << id << " has no material named " << name); AKANTU_DEBUG_OUT(); return *materials[it->second]; } /* -------------------------------------------------------------------------- */ inline void SolidMechanicsModel::setMaterialSelector(MaterialSelector & selector) { if (is_default_material_selector) delete material_selector; material_selector = &selector; is_default_material_selector = false; } -/* -------------------------------------------------------------------------- */ -inline FEEngine & SolidMechanicsModel::getFEEngineBoundary(const ID & name) { - return dynamic_cast<FEEngine &>( - getFEEngineClassBoundary<MyFEEngineType>(name)); -} /* -------------------------------------------------------------------------- */ inline void SolidMechanicsModel::splitElementByMaterial( const Array<Element> & elements, Array<Element> * elements_per_mat) const { ElementType current_element_type = _not_defined; GhostType current_ghost_type = _casper; const Array<UInt> * mat_indexes = NULL; const Array<UInt> * mat_loc_num = NULL; Array<Element>::const_iterator<Element> it = elements.begin(); Array<Element>::const_iterator<Element> end = elements.end(); for (; it != end; ++it) { Element el = *it; if (el.type != current_element_type || el.ghost_type != current_ghost_type) { current_element_type = el.type; current_ghost_type = el.ghost_type; mat_indexes = &(this->material_index(el.type, el.ghost_type)); mat_loc_num = &(this->material_local_numbering(el.type, el.ghost_type)); } UInt old_id = el.element; el.element = (*mat_loc_num)(old_id); elements_per_mat[(*mat_indexes)(old_id)].push_back(el); } } /* -------------------------------------------------------------------------- */ inline UInt SolidMechanicsModel::getNbData(const Array<Element> & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes_per_element = 0; Array<Element>::const_iterator<Element> it = elements.begin(); Array<Element>::const_iterator<Element> end = elements.end(); for (; it != end; ++it) { const Element & el = *it; nb_nodes_per_element += Mesh::getNbNodesPerElement(el.type); } switch (tag) { case _gst_material_id: { size += elements.getSize() * sizeof(UInt); break; } case _gst_smm_mass: { size += nb_nodes_per_element * sizeof(Real) * Model::spatial_dimension; // mass vector break; } case _gst_smm_for_gradu: { size += nb_nodes_per_element * Model::spatial_dimension * sizeof(Real); // displacement break; } case _gst_smm_boundary: { // force, displacement, boundary size += nb_nodes_per_element * Model::spatial_dimension * (2 * sizeof(Real) + sizeof(bool)); break; } case _gst_for_dump: { // displacement, velocity, acceleration, residual, force size += nb_nodes_per_element * Model::spatial_dimension * sizeof(Real) * 5; break; } default: {} } if (tag != _gst_material_id) { Array<Element> * elements_per_mat = new Array<Element>[materials.size()]; this->splitElementByMaterial(elements, elements_per_mat); for (UInt i = 0; i < materials.size(); ++i) { size += materials[i]->getNbData(elements_per_mat[i], tag); } delete[] elements_per_mat; } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ inline void SolidMechanicsModel::packData(CommunicationBuffer & buffer, const Array<Element> & elements, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); switch (tag) { case _gst_material_id: { this->packElementalDataHelper(material_index, buffer, elements, false, getFEEngine()); break; } case _gst_smm_mass: { packNodalDataHelper(*mass, buffer, elements, mesh); break; } case _gst_smm_for_gradu: { packNodalDataHelper(*displacement, buffer, elements, mesh); break; } case _gst_for_dump: { packNodalDataHelper(*displacement, buffer, elements, mesh); packNodalDataHelper(*velocity, buffer, elements, mesh); packNodalDataHelper(*acceleration, buffer, elements, mesh); packNodalDataHelper(*internal_force, buffer, elements, mesh); packNodalDataHelper(*external_force, buffer, elements, mesh); break; } case _gst_smm_boundary: { packNodalDataHelper(*external_force, buffer, elements, mesh); packNodalDataHelper(*velocity, buffer, elements, mesh); packNodalDataHelper(*blocked_dofs, buffer, elements, mesh); break; } default: {} } if (tag != _gst_material_id) { Array<Element> * elements_per_mat = new Array<Element>[materials.size()]; splitElementByMaterial(elements, elements_per_mat); for (UInt i = 0; i < materials.size(); ++i) { materials[i]->packData(buffer, elements_per_mat[i], tag); } delete[] elements_per_mat; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer, const Array<Element> & elements, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); switch (tag) { case _gst_material_id: { unpackElementalDataHelper(material_index, buffer, elements, false, getFEEngine()); break; } case _gst_smm_mass: { unpackNodalDataHelper(*mass, buffer, elements, mesh); break; } case _gst_smm_for_gradu: { unpackNodalDataHelper(*displacement, buffer, elements, mesh); break; } case _gst_for_dump: { unpackNodalDataHelper(*displacement, buffer, elements, mesh); unpackNodalDataHelper(*velocity, buffer, elements, mesh); unpackNodalDataHelper(*acceleration, buffer, elements, mesh); unpackNodalDataHelper(*internal_force, buffer, elements, mesh); unpackNodalDataHelper(*external_force, buffer, elements, mesh); break; } case _gst_smm_boundary: { unpackNodalDataHelper(*external_force, buffer, elements, mesh); unpackNodalDataHelper(*velocity, buffer, elements, mesh); unpackNodalDataHelper(*blocked_dofs, buffer, elements, mesh); break; } default: {} } if (tag != _gst_material_id) { Array<Element> * elements_per_mat = new Array<Element>[materials.size()]; splitElementByMaterial(elements, elements_per_mat); for (UInt i = 0; i < materials.size(); ++i) { materials[i]->unpackData(buffer, elements_per_mat[i], tag); } delete[] elements_per_mat; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline UInt SolidMechanicsModel::getNbData(const Array<UInt> & dofs, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); UInt size = 0; // UInt nb_nodes = mesh.getNbNodes(); switch (tag) { case _gst_smm_uv: { size += sizeof(Real) * Model::spatial_dimension * 2; break; } case _gst_smm_res: { size += sizeof(Real) * Model::spatial_dimension; break; } case _gst_smm_mass: { size += sizeof(Real) * Model::spatial_dimension; break; } case _gst_for_dump: { size += sizeof(Real) * Model::spatial_dimension * 5; break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size * dofs.getSize(); } /* -------------------------------------------------------------------------- */ inline void SolidMechanicsModel::packData(CommunicationBuffer & buffer, const Array<UInt> & dofs, const SynchronizationTag & tag) const { AKANTU_DEBUG_IN(); switch (tag) { case _gst_smm_uv: { packDOFDataHelper(*displacement, buffer, dofs); packDOFDataHelper(*velocity, buffer, dofs); break; } case _gst_smm_res: { packDOFDataHelper(*internal_force, buffer, dofs); break; } case _gst_smm_mass: { packDOFDataHelper(*mass, buffer, dofs); break; } case _gst_for_dump: { packDOFDataHelper(*displacement, buffer, dofs); packDOFDataHelper(*velocity, buffer, dofs); packDOFDataHelper(*acceleration, buffer, dofs); packDOFDataHelper(*internal_force, buffer, dofs); packDOFDataHelper(*external_force, buffer, dofs); break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline void SolidMechanicsModel::unpackData(CommunicationBuffer & buffer, const Array<UInt> & dofs, const SynchronizationTag & tag) { AKANTU_DEBUG_IN(); switch (tag) { case _gst_smm_uv: { unpackDOFDataHelper(*displacement, buffer, dofs); unpackDOFDataHelper(*velocity, buffer, dofs); break; } case _gst_smm_res: { unpackDOFDataHelper(*internal_force, buffer, dofs); break; } case _gst_smm_mass: { unpackDOFDataHelper(*mass, buffer, dofs); break; } case _gst_for_dump: { unpackDOFDataHelper(*displacement, buffer, dofs); unpackDOFDataHelper(*velocity, buffer, dofs); unpackDOFDataHelper(*acceleration, buffer, dofs); unpackDOFDataHelper(*internal_force, buffer, dofs); unpackDOFDataHelper(*external_force, buffer, dofs); break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ } // namespace akantu #endif /* __AKANTU_SOLID_MECHANICS_MODEL_INLINE_IMPL_CC__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_mass.cc b/src/model/solid_mechanics/solid_mechanics_model_mass.cc index ec921354c..fa5fcca40 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_mass.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_mass.cc @@ -1,170 +1,172 @@ /** * @file solid_mechanics_model_mass.cc * * @author Nicolas Richart <nicolas.richart@epfl.ch> * * @date creation: Tue Oct 05 2010 * @date last modification: Fri Oct 16 2015 * * @brief function handling mass computation * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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 "integrator_gauss.hh" #include "material.hh" #include "model_solver.hh" +#include "shape_lagrange.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ namespace akantu { class ComputeRhoFunctor { public: explicit ComputeRhoFunctor(const SolidMechanicsModel & model) : model(model){}; void operator()(Matrix<Real> & rho, const Element & element) const { const Array<UInt> & mat_indexes = model.getMaterialByElement(element.type, element.ghost_type); Real mat_rho = model.getMaterial(mat_indexes(element.element)).getParam("rho"); rho.set(mat_rho); } private: const SolidMechanicsModel & model; }; /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMassLumped() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); if (this->mass == NULL) { std::stringstream sstr_mass; sstr_mass << id << ":mass"; mass = &(alloc<Real>(sstr_mass.str(), nb_nodes, Model::spatial_dimension, 0)); } else { mass->clear(); } if (!this->getDOFManager().hasLumpedMatrix("M")) { this->getDOFManager().getNewLumpedMatrix("M"); } this->getDOFManager().clearLumpedMatrix("M"); assembleMassLumped(_not_ghost); assembleMassLumped(_ghost); this->getDOFManager().getLumpedMatrixPerDOFs("displacement", "M", *(this->mass)); /// for not connected nodes put mass to one in order to avoid #if !defined(AKANTU_NDEBUG) bool has_unconnected_nodes = false; auto mass_it = mass->begin_reinterpret(mass->getSize() * mass->getNbComponent()); auto mass_end = mass->end_reinterpret(mass->getSize() * mass->getNbComponent()); for (; mass_it != mass_end; ++mass_it) { if (std::abs(*mass_it) < std::numeric_limits<Real>::epsilon() || Math::isnan(*mass_it)) { has_unconnected_nodes = true; break; } } if (has_unconnected_nodes) AKANTU_DEBUG_WARNING("There are nodes that seem to not be connected to any " "elements, beware that they have lumped mass of 0."); #endif this->synchronize(_gst_smm_mass); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMass() { AKANTU_DEBUG_IN(); if (!this->getDOFManager().hasMatrix("M")) { this->getDOFManager().getNewMatrix("M", "J"); } this->getDOFManager().clearMatrix("M"); assembleMass(_not_ghost); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMassLumped(GhostType ghost_type) { AKANTU_DEBUG_IN(); MyFEEngineType & fem = getFEEngineClass<MyFEEngineType>(); ComputeRhoFunctor compute_rho(*this); for (auto type : mesh.elementTypes(Model::spatial_dimension, ghost_type)) { fem.assembleFieldLumped(compute_rho, "M", "displacement", this->getDOFManager(), type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMass(GhostType ghost_type) { AKANTU_DEBUG_IN(); MyFEEngineType & fem = getFEEngineClass<MyFEEngineType>(); ComputeRhoFunctor compute_rho(*this); for (auto type : mesh.elementTypes(Model::spatial_dimension, ghost_type)) { fem.assembleFieldMatrix(compute_rho, "M", "displacement", this->getDOFManager(), type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMassLumped(const Array<UInt> &) { AKANTU_DEBUG_IN(); assembleMassLumped(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMass(const Array<UInt> &) { AKANTU_DEBUG_IN(); assembleMass(); AKANTU_DEBUG_OUT(); } } // namespace akantu diff --git a/src/synchronizer/master_element_info_per_processor.cc b/src/synchronizer/master_element_info_per_processor.cc index bc6a866ca..15742b202 100644 --- a/src/synchronizer/master_element_info_per_processor.cc +++ b/src/synchronizer/master_element_info_per_processor.cc @@ -1,485 +1,482 @@ /** * @file master_element_info_per_processor.cc * * @author Nicolas Richart <nicolas.richart@epfl.ch> * * @date Fri Mar 11 14:57:13 2016 * * @brief * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see <http://www.gnu.org/licenses/>. * */ /* -------------------------------------------------------------------------- */ #include "aka_iterators.hh" +#include "mesh_iterators.hh" #include "element_group.hh" #include "element_info_per_processor.hh" #include "element_synchronizer.hh" #include "mesh_utils.hh" #include "static_communicator.hh" /* -------------------------------------------------------------------------- */ #include <algorithm> #include <iostream> #include <map> #include <tuple> /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ MasterElementInfoPerProc::MasterElementInfoPerProc( ElementSynchronizer & synchronizer, UInt message_cnt, UInt root, ElementType type, const MeshPartition & partition) : ElementInfoPerProc(synchronizer, message_cnt, root, type), partition(partition), all_nb_local_element(nb_proc, 0), all_nb_ghost_element(nb_proc, 0), all_nb_element_to_send(nb_proc, 0) { Vector<UInt> size(5); size(0) = (UInt)type; if (type != _not_defined) { nb_nodes_per_element = Mesh::getNbNodesPerElement(type); nb_element = mesh.getNbElement(type); const Array<UInt> & partition_num = this->partition.getPartition(this->type, _not_ghost); const CSR<UInt> & ghost_partition = this->partition.getGhostPartitionCSR()(this->type, _not_ghost); for (UInt el = 0; el < nb_element; ++el) { this->all_nb_local_element[partition_num(el)]++; for (CSR<UInt>::const_iterator part = ghost_partition.begin(el); part != ghost_partition.end(el); ++part) { this->all_nb_ghost_element[*part]++; } this->all_nb_element_to_send[partition_num(el)] += ghost_partition.getNbCols(el) + 1; } /// tag info std::vector<std::string> tag_names; this->getMeshData().getTagNames(tag_names, type); this->nb_tags = tag_names.size(); size(4) = nb_tags; for (UInt p = 0; p < nb_proc; ++p) { if (p != root) { size(1) = this->all_nb_local_element[p]; size(2) = this->all_nb_ghost_element[p]; size(3) = this->all_nb_element_to_send[p]; AKANTU_DEBUG_INFO( "Sending connectivities informations to proc " << p << " TAG(" << Tag::genTag(this->rank, this->message_count, Tag::_SIZES) << ")"); comm.send(size, p, Tag::genTag(this->rank, this->message_count, Tag::_SIZES)); } else { this->nb_local_element = this->all_nb_local_element[p]; this->nb_ghost_element = this->all_nb_ghost_element[p]; } } } else { for (UInt p = 0; p < this->nb_proc; ++p) { if (p != this->root) { AKANTU_DEBUG_INFO( "Sending empty connectivities informations to proc " << p << " TAG(" << Tag::genTag(this->rank, this->message_count, Tag::_SIZES) << ")"); comm.send(size, p, Tag::genTag(this->rank, this->message_count, Tag::_SIZES)); } } } } /* ------------------------------------------------------------------------ */ void MasterElementInfoPerProc::synchronizeConnectivities() { const Array<UInt> & partition_num = this->partition.getPartition(this->type, _not_ghost); const CSR<UInt> & ghost_partition = this->partition.getGhostPartitionCSR()(this->type, _not_ghost); std::vector<Array<UInt>> buffers(this->nb_proc); auto conn_it = this->mesh.getConnectivity(this->type, _not_ghost) .begin(this->nb_nodes_per_element); auto conn_end = this->mesh.getConnectivity(this->type, _not_ghost) .end(this->nb_nodes_per_element); /// copying the local connectivity auto part_it = partition_num.begin(); for (; conn_it != conn_end; ++conn_it, ++part_it) { const auto & conn = *conn_it; for (UInt i = 0; i < conn.size(); ++i) { buffers[*part_it].push_back(conn[i]); } } /// copying the connectivity of ghost element conn_it = this->mesh.getConnectivity(this->type, _not_ghost) .begin(this->nb_nodes_per_element); for (UInt el = 0; conn_it != conn_end; ++el, ++conn_it) { for (auto part = ghost_partition.begin(el); part != ghost_partition.end(el); ++part) { UInt proc = *part; const Vector<UInt> & conn = *conn_it; for (UInt i = 0; i < conn.size(); ++i) { buffers[proc].push_back(conn[i]); } } } #ifndef AKANTU_NDEBUG for (UInt p = 0; p < this->nb_proc; ++p) { UInt size = this->nb_nodes_per_element * (this->all_nb_local_element[p] + this->all_nb_ghost_element[p]); AKANTU_DEBUG_ASSERT( buffers[p].getSize() == size, "The connectivity data packed in the buffer are not correct"); } #endif /// send all connectivity and ghost information to all processors std::vector<CommunicationRequest> requests; for (UInt p = 0; p < this->nb_proc; ++p) { if (p != this->root) { AKANTU_DEBUG_INFO( "Sending connectivities to proc " << p << " TAG(" << Tag::genTag(this->rank, this->message_count, Tag::_CONNECTIVITY) << ")"); requests.push_back(comm.asyncSend( buffers[p], p, Tag::genTag(this->rank, this->message_count, Tag::_CONNECTIVITY))); } } Array<UInt> & old_nodes = this->getNodesGlobalIds(); /// create the renumbered connectivity AKANTU_DEBUG_INFO("Renumbering local connectivities"); MeshUtils::renumberMeshNodes(mesh, buffers[root], all_nb_local_element[root], all_nb_ghost_element[root], type, old_nodes); comm.waitAll(requests); comm.freeCommunicationRequest(requests); } /* ------------------------------------------------------------------------ */ void MasterElementInfoPerProc::synchronizePartitions() { const Array<UInt> & partition_num = this->partition.getPartition(this->type, _not_ghost); const CSR<UInt> & ghost_partition = this->partition.getGhostPartitionCSR()(this->type, _not_ghost); std::vector<Array<UInt>> buffers(this->partition.getNbPartition()); /// splitting the partition information to send them to processors Vector<UInt> count_by_proc(nb_proc, 0); for (UInt el = 0; el < nb_element; ++el) { UInt proc = partition_num(el); buffers[proc].push_back(ghost_partition.getNbCols(el)); UInt i(0); for (CSR<UInt>::const_iterator part = ghost_partition.begin(el); part != ghost_partition.end(el); ++part, ++i) { buffers[proc].push_back(*part); } } for (UInt el = 0; el < nb_element; ++el) { UInt i(0); for (CSR<UInt>::const_iterator part = ghost_partition.begin(el); part != ghost_partition.end(el); ++part, ++i) { buffers[*part].push_back(partition_num(el)); } } #ifndef AKANTU_NDEBUG for (UInt p = 0; p < this->nb_proc; ++p) { AKANTU_DEBUG_ASSERT( buffers[p].getSize() == (this->all_nb_ghost_element[p] + this->all_nb_element_to_send[p]), "Data stored in the buffer are most probably wrong"); } #endif std::vector<CommunicationRequest> requests; /// last data to compute the communication scheme for (UInt p = 0; p < this->nb_proc; ++p) { if (p != this->root) { AKANTU_DEBUG_INFO( "Sending partition informations to proc " << p << " TAG(" << Tag::genTag(this->rank, this->message_count, Tag::_PARTITIONS) << ")"); requests.push_back(comm.asyncSend( buffers[p], p, Tag::genTag(this->rank, this->message_count, Tag::_PARTITIONS))); } } if (Mesh::getSpatialDimension(this->type) == this->mesh.getSpatialDimension()) { AKANTU_DEBUG_INFO("Creating communications scheme"); this->fillCommunicationScheme(buffers[this->rank]); } comm.waitAll(requests); comm.freeCommunicationRequest(requests); } /* -------------------------------------------------------------------------- */ void MasterElementInfoPerProc::synchronizeTags() { AKANTU_DEBUG_IN(); if (this->nb_tags == 0) { AKANTU_DEBUG_OUT(); return; } UInt mesh_data_sizes_buffer_length; MeshData & mesh_data = this->getMeshData(); /// tag info std::vector<std::string> tag_names; mesh_data.getTagNames(tag_names, type); // Make sure the tags are sorted (or at least not in random order), // because they come from a map !! std::sort(tag_names.begin(), tag_names.end()); // Sending information about the tags in mesh_data: name, data type and // number of components of the underlying array associated to the current // type DynamicCommunicationBuffer mesh_data_sizes_buffer; std::vector<std::string>::const_iterator names_it = tag_names.begin(); std::vector<std::string>::const_iterator names_end = tag_names.end(); for (; names_it != names_end; ++names_it) { mesh_data_sizes_buffer << *names_it; mesh_data_sizes_buffer << mesh_data.getTypeCode(*names_it); mesh_data_sizes_buffer << mesh_data.getNbComponent(*names_it, type); } mesh_data_sizes_buffer_length = mesh_data_sizes_buffer.getSize(); AKANTU_DEBUG_INFO( "Broadcasting the size of the information about the mesh data tags: (" << mesh_data_sizes_buffer_length << ")."); comm.broadcast(mesh_data_sizes_buffer_length, root); AKANTU_DEBUG_INFO( "Broadcasting the information about the mesh data tags, addr " << (void *)mesh_data_sizes_buffer.storage()); if (mesh_data_sizes_buffer_length != 0) comm.broadcast(mesh_data_sizes_buffer, root); if (mesh_data_sizes_buffer_length != 0) { // Sending the actual data to each processor DynamicCommunicationBuffer * buffers = new DynamicCommunicationBuffer[nb_proc]; std::vector<std::string>::const_iterator names_it = tag_names.begin(); std::vector<std::string>::const_iterator names_end = tag_names.end(); // Loop over each tag for the current type for (; names_it != names_end; ++names_it) { // Type code of the current tag (i.e. the tag named *names_it) this->fillTagBuffer(buffers, *names_it); } std::vector<CommunicationRequest> requests; for (UInt p = 0; p < nb_proc; ++p) { if (p != root) { AKANTU_DEBUG_INFO("Sending " << buffers[p].getSize() << " bytes of mesh data to proc " << p << " TAG(" << Tag::genTag(this->rank, this->message_count, Tag::_MESH_DATA) << ")"); requests.push_back(comm.asyncSend( buffers[p], p, Tag::genTag(this->rank, this->message_count, Tag::_MESH_DATA))); } } names_it = tag_names.begin(); // Loop over each tag for the current type for (; names_it != names_end; ++names_it) { // Reinitializing the mesh data on the master this->fillMeshData(buffers[root], *names_it, mesh_data.getTypeCode(*names_it), mesh_data.getNbComponent(*names_it, type)); } comm.waitAll(requests); comm.freeCommunicationRequest(requests); requests.clear(); delete[] buffers; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template <typename T> void MasterElementInfoPerProc::fillTagBufferTemplated( DynamicCommunicationBuffer * buffers, const std::string & tag_name) { MeshData & mesh_data = this->getMeshData(); const Array<T> & data = mesh_data.getElementalDataArray<T>(tag_name, type); const Array<UInt> & partition_num = this->partition.getPartition(this->type, _not_ghost); const CSR<UInt> & ghost_partition = this->partition.getGhostPartitionCSR()(this->type, _not_ghost); // Not possible to use the iterator because it potentially triggers the // creation of complex // type templates (such as akantu::Vector< std::vector<Element> > which don't // implement the right interface // (e.g. operator<< in that case). // typename Array<T>::template const_iterator< Vector<T> > data_it = // data.begin(data.getNbComponent()); // typename Array<T>::template const_iterator< Vector<T> > data_end = // data.end(data.getNbComponent()); const T * data_it = data.storage(); const T * data_end = data.storage() + data.getSize() * data.getNbComponent(); const UInt * part = partition_num.storage(); /// copying the data, element by element for (; data_it != data_end; ++part) { for (UInt j(0); j < data.getNbComponent(); ++j, ++data_it) { buffers[*part] << *data_it; } } data_it = data.storage(); /// copying the data for the ghost element for (UInt el(0); data_it != data_end; data_it += data.getNbComponent(), ++el) { CSR<UInt>::const_iterator it = ghost_partition.begin(el); CSR<UInt>::const_iterator end = ghost_partition.end(el); for (; it != end; ++it) { UInt proc = *it; for (UInt j(0); j < data.getNbComponent(); ++j) { buffers[proc] << data_it[j]; } } } } /* -------------------------------------------------------------------------- */ void MasterElementInfoPerProc::fillTagBuffer( DynamicCommunicationBuffer * buffers, const std::string & tag_name) { MeshData & mesh_data = this->getMeshData(); #define AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA(r, extra_param, elem) \ case BOOST_PP_TUPLE_ELEM(2, 0, elem): { \ this->fillTagBufferTemplated<BOOST_PP_TUPLE_ELEM(2, 1, elem)>(buffers, \ tag_name); \ break; \ } MeshDataTypeCode data_type_code = mesh_data.getTypeCode(tag_name); switch (data_type_code) { BOOST_PP_SEQ_FOR_EACH(AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA, , AKANTU_MESH_DATA_TYPES) default: AKANTU_DEBUG_ERROR("Could not obtain the type of tag" << tag_name << "!"); break; } #undef AKANTU_DISTRIBUTED_SYNHRONIZER_TAG_DATA } /* -------------------------------------------------------------------------- */ void MasterElementInfoPerProc::synchronizeGroups() { AKANTU_DEBUG_IN(); DynamicCommunicationBuffer * buffers = new DynamicCommunicationBuffer[nb_proc]; using ElementToGroup = std::vector<std::vector<std::string>>; ElementToGroup element_to_group; element_to_group.resize(nb_element); - auto egi = mesh.element_group_begin(); - auto ege = mesh.element_group_end(); - for (; egi != ege; ++egi) { - ElementGroup & eg = *(egi->second); - - std::string name = egi->first; + for (auto & eg : ElementGroupsIterable(mesh)) { + const std::string & name = eg.getName(); for (const auto & element : eg.getElements(type, _not_ghost)) { element_to_group[element].push_back(name); } auto eit = eg.begin(type, _not_ghost); if (eit != eg.end(type, _not_ghost)) const_cast<Array<UInt> &>(eg.getElements(type)).empty(); } const auto & partition_num = this->partition.getPartition(this->type, _not_ghost); const auto & ghost_partition = this->partition.getGhostPartitionCSR()(this->type, _not_ghost); /// copying the data, element by element ElementToGroup::const_iterator data_it = element_to_group.begin(); ElementToGroup::const_iterator data_end = element_to_group.end(); for (auto pair : zip(partition_num, element_to_group)) { buffers[std::get<0>(pair)] << std::get<1>(pair); } data_it = element_to_group.begin(); /// copying the data for the ghost element for (UInt el(0); data_it != data_end; ++data_it, ++el) { CSR<UInt>::const_iterator it = ghost_partition.begin(el); CSR<UInt>::const_iterator end = ghost_partition.end(el); for (; it != end; ++it) { UInt proc = *it; buffers[proc] << *data_it; } } std::vector<CommunicationRequest> requests; for (UInt p = 0; p < this->nb_proc; ++p) { if (p == this->rank) continue; AKANTU_DEBUG_INFO("Sending element groups to proc " << p << " TAG(" << Tag::genTag(this->rank, p, Tag::_ELEMENT_GROUP) << ")"); requests.push_back(comm.asyncSend( buffers[p], p, Tag::genTag(this->rank, p, Tag::_ELEMENT_GROUP))); } this->fillElementGroupsFromBuffer(buffers[this->rank]); comm.waitAll(requests); comm.freeCommunicationRequest(requests); requests.clear(); delete[] buffers; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ } // namespace akantu diff --git a/test/test_common/test_arange_iterator.cc b/test/test_common/test_arange_iterator.cc index 2228067a6..94463e101 100644 --- a/test/test_common/test_arange_iterator.cc +++ b/test/test_common/test_arange_iterator.cc @@ -1,48 +1,48 @@ /** * @file test_arange_iterator.cc * * @author Nicolas Richart * * @date creation Fri Aug 11 2017 * * @brief A Documented file. * * @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_iterators.hh" /* -------------------------------------------------------------------------- */ #include <iostream> /* -------------------------------------------------------------------------- */ using namespace akantu; int main() { for (auto i : arange(10)) std::cout << i << std::endl; - for (auto i : arange(-1, 10, 10)) + for (auto i : arange(1, 22, 2)) std::cout << i << std::endl; - for (auto i : arange(-1, -10, -1)) - std::cout << i << std::endl; + for (auto && i : zip(arange(-1, -10, -1), arange(1, 18, 2))) + std::cout << std::get<0>(i) << " " << std::get<1>(i) << std::endl; return 0; } diff --git a/test/test_fe_engine/test_inverse_map.cc b/test/test_fe_engine/test_inverse_map.cc index 9ecdcb234..daeac8134 100644 --- a/test/test_fe_engine/test_inverse_map.cc +++ b/test/test_fe_engine/test_inverse_map.cc @@ -1,106 +1,110 @@ /** * @file test_inverse_map.cc * * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch> * * @date creation: Fri Sep 03 2010 * @date last modification: Thu Oct 15 2015 * * @brief test of the fem class * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see <http://www.gnu.org/licenses/>. * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" +#include "aka_iterators.hh" #include "fe_engine.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" /* -------------------------------------------------------------------------- */ -#include <cstdlib> #include <iostream> /* -------------------------------------------------------------------------- */ using namespace akantu; int main(int argc, char * argv[]) { akantu::initialize(argc, argv); - debug::setDebugLevel(dblTest); - const ElementType type = TYPE; - UInt dim = ElementClass<type>::getSpatialDimension(); + const auto type = TYPE; + auto dim = ElementClass<type>::getSpatialDimension(); Mesh my_mesh(dim); - - const Vector<Real> & lower = my_mesh.getLowerBounds(); - const Vector<Real> & upper = my_mesh.getUpperBounds(); - std::stringstream meshfilename; meshfilename << type << ".msh"; + my_mesh.read(meshfilename.str()); + const auto & lower = my_mesh.getLowerBounds(); + const auto & upper = my_mesh.getUpperBounds(); UInt nb_elements = my_mesh.getNbElement(type); - /// - FEEngineTemplate<IntegratorGauss, ShapeLagrange> * fem = - new FEEngineTemplate<IntegratorGauss, ShapeLagrange>(my_mesh, dim, - "my_fem"); - fem->initShapeFunctions(); + auto fem = std::make_unique<FEEngineTemplate<IntegratorGauss, ShapeLagrange>>( + my_mesh, dim, "my_fem"); - UInt nb_quad_points = fem->getNbIntegrationPoints(type); + fem->initShapeFunctions(); + Matrix<Real> quad = GaussIntegrationElement<type>::getQuadraturePoints(); /// get the quadrature points coordinates - Array<Real> coord_on_quad(nb_quad_points * nb_elements, + Array<Real> coord_on_quad(quad.cols() * nb_elements, my_mesh.getSpatialDimension(), "coord_on_quad"); fem->interpolateOnIntegrationPoints(my_mesh.getNodes(), coord_on_quad, my_mesh.getSpatialDimension(), type); - /// loop over the quadrature points - Array<Real>::iterator<Vector<Real>> it = coord_on_quad.begin(dim); Vector<Real> natural_coords(dim); - Matrix<Real> quad = GaussIntegrationElement<type>::getQuadraturePoints(); - - for (UInt el = 0; el < nb_elements; ++el) { - for (UInt q = 0; q < nb_quad_points; ++q) { - fem->inverseMap(*it, el, type, natural_coords); - for (UInt i = 0; i < dim; ++i) { - __attribute__((unused)) const Real eps = 1e-13; - AKANTU_DEBUG_ASSERT( - std::abs((natural_coords(i) - quad(i, q)) / (upper(i) - lower(i))) < - eps, - "real coordinates inversion test failed:" - << natural_coords(i) << " - " << quad(i, q) << " = " - << (natural_coords(i) - quad(i, q)) / (upper(i) - lower(i))); + /// loop over the quadrature points + auto it = coord_on_quad.begin_reinterpret(dim, quad.cols(), nb_elements); + auto end = coord_on_quad.end_reinterpret(dim, quad.cols(), nb_elements); + + auto length = (upper - lower).norm<L_inf>(); + + auto eps = 5e-12; + UInt el = 0; + for (; it != end; ++it, ++el) { + const auto & quads_coords = *it; + for (auto q : arange(quad.cols())) { + Vector<Real> quad_coord = quads_coords(q); + Vector<Real> ref_quad_coord = quad(q); + fem->inverseMap(quad_coord, el, type, natural_coords); + + auto dis_normalized = ref_quad_coord.distance(natural_coords) / length; + if (not(dis_normalized < eps)) { + std::cout << "Real coordinates inversion test failed : " + << std::scientific << natural_coords << " - " + << ref_quad_coord << " / " << length << " = " << dis_normalized << std::endl; + return 1; } - ++it; + + // std::cout << "Real coordinates inversion : " << std::scientific + // << natural_coords << " = " << ref_quad_coord << " (" + // << dis_normalized << ")" << std::endl; } } std::cout << "inverse completed over " << nb_elements << " elements" << std::endl; - delete fem; finalize(); - return EXIT_SUCCESS; + return 0; } diff --git a/test/test_fe_engine/test_mesh_boundary.cc b/test/test_fe_engine/test_mesh_boundary.cc index 2a0c0d042..9e34b4e99 100644 --- a/test/test_fe_engine/test_mesh_boundary.cc +++ b/test/test_fe_engine/test_mesh_boundary.cc @@ -1,78 +1,64 @@ /** * @file test_mesh_boundary.cc * * @author Dana Christen <dana.christen@gmail.com> * * @date creation: Fri May 03 2013 * @date last modification: Mon Jul 13 2015 * * @brief Thest the element groups * * @section LICENSE * * Copyright (©) 2014, 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see <http://www.gnu.org/licenses/>. * */ - -#include <iostream> -#include <sstream> -#include "aka_common.hh" +/* -------------------------------------------------------------------------- */ #include "mesh.hh" +/* -------------------------------------------------------------------------- */ +#include <iostream> +/* -------------------------------------------------------------------------- */ using namespace akantu; -/* -------------------------------------------------------------------------- */ - -int main(int argc, char* argv[]) { +int main(int argc, char * argv[]) { UInt spatialDimension(3); akantu::initialize(argc, argv); Mesh mesh(spatialDimension, "mesh_names"); std::cout << "Loading the mesh." << std::endl; - // mesh.read("./cube_physical_names.msh"); mesh.read("./cube_physical_names.msh"); - std::stringstream sstr; std::cout << "Examining mesh:" << std::endl; // Inspection of the number of boundaries - __attribute__ ((unused)) UInt nb_boundaries= mesh.getNbElementGroups(); - AKANTU_DEBUG_INFO(nb_boundaries << " boundaries advertised initially by Mesh."); - - AKANTU_DEBUG_INFO("Building boundaries"); - - // Two methods: either building using data loaded from the mesh file in MeshData - // or build with automatic numbering - mesh.createGroupsFromMeshData<std::string>("physical_names"); - - // Second inspection of the number of boundaries (should not be 0) - nb_boundaries = mesh.getNbElementGroups(); - + UInt nb_boundaries = mesh.getNbElementGroups(spatialDimension - 1); AKANTU_DEBUG_INFO(nb_boundaries << " boundaries advertised by Mesh."); - AKANTU_DEBUG_ASSERT(nb_boundaries != 0, "No boundary detected!"); + if (nb_boundaries == 0) { + std::cout << "No boundary detected!" << std::endl; + return 1; + } - std::cout << (*dynamic_cast<GroupManager*>(&mesh)) << std::endl; + std::cout << (*dynamic_cast<GroupManager *>(&mesh)) << std::endl; akantu::finalize(); - return EXIT_SUCCESS; + return 0; } - - diff --git a/test/test_mesh_utils/test_mesh_iterators.cc b/test/test_mesh_utils/test_mesh_iterators.cc index d9a01e714..111602919 100644 --- a/test/test_mesh_utils/test_mesh_iterators.cc +++ b/test/test_mesh_utils/test_mesh_iterators.cc @@ -1,56 +1,67 @@ /** * @file test_mesh_iterators.cc * * @author Nicolas Richart * * @date creation Wed Aug 09 2017 * * @brief Test the mesh iterators * * @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_iterators.hh" #include "element_group.hh" #include "mesh.hh" #include "mesh_iterators.hh" #include "node_group.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; int main(int argc, char * argv[]) { initialize(argc, argv); Mesh mesh(3); const Mesh & cmesh = mesh; mesh.read("iterators_mesh.msh"); - for (auto && element_group : MeshElementGroups(mesh)) { + for (auto && element_group : ElementGroupsIterable(mesh)) { std::cout << element_group.getName() << std::endl; } - for (auto && node_group : MeshNodeGroups(cmesh)) { + for (auto && node_group : NodeGroupsIterable(cmesh)) { std::cout << node_group.getName() << std::endl; } + for (auto && element_group : counting(ElementGroupsIterable(mesh))) { + std::cout << std::get<0>(element_group) << " " << std::get<1>(element_group).getName() << std::endl; + + } + + // for (auto && node_group : + // counting(NodeGroupsIterable(cmesh))) { + // std::cout << std::get<0>(node_group) << " " << std::get<1>(node_group).getName() << std::endl; + // } + finalize(); return EXIT_SUCCESS; } diff --git a/test/test_mesh_utils/test_purify_mesh.cc b/test/test_mesh_utils/test_purify_mesh.cc index 72de7434d..a9c19be5a 100644 --- a/test/test_mesh_utils/test_purify_mesh.cc +++ b/test/test_mesh_utils/test_purify_mesh.cc @@ -1,63 +1,62 @@ /** * @file test_purify_mesh.cc * * @author Nicolas Richart <nicolas.richart@epfl.ch> * * @date creation: Thu Jul 15 2010 * @date last modification: Sun Oct 19 2014 * * @brief Test the purifyMesh function from MeshUtils * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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 "mesh.hh" #include "mesh_utils.hh" #include "mesh_io.hh" using namespace akantu; int main(int argc, char *argv[]) { akantu::initialize(argc, argv); Mesh mesh(2); - MeshIOMSH mesh_io; - mesh_io.read("purify_mesh.msh", mesh); + mesh.read("purify_mesh.msh"); MeshUtils::purifyMesh(mesh); - mesh_io.write("purify_mesh_after.msh", mesh); + mesh.write("purify_mesh_after.msh"); if(mesh.getNbNodes() != 21) AKANTU_DEBUG_ERROR("The purified mesh does not contain the good number of nodes."); if(mesh.getNbElement(_quadrangle_8) != 4) AKANTU_DEBUG_ERROR("The purified mesh does not contain the good number of element."); akantu::finalize(); return EXIT_SUCCESS; } diff --git a/test/test_model/test_solid_mechanics_model/test_materials/test_interpolate_stress.cc b/test/test_model/test_solid_mechanics_model/test_materials/test_interpolate_stress.cc index 952d1bb85..264dfe49a 100644 --- a/test/test_model/test_solid_mechanics_model/test_materials/test_interpolate_stress.cc +++ b/test/test_model/test_solid_mechanics_model/test_materials/test_interpolate_stress.cc @@ -1,181 +1,182 @@ /** * @file test_interpolate_stress.cc * * @author Marco Vocialta <marco.vocialta@epfl.ch> * * @date creation: Thu Jun 07 2012 * @date last modification: Thu Oct 15 2015 * * @brief Test for the stress interpolation function * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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 <fstream> #include <iostream> #include <limits> /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" - +#include "integrator_gauss.hh" +#include "shape_lagrange.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; Real function(Real x, Real y, Real z) { return 100. + 2. * x + 3. * y + 4 * z; } int main(int argc, char * argv[]) { initialize("../material.dat", argc, argv); debug::setDebugLevel(dblWarning); const UInt spatial_dimension = 3; const ElementType type = _tetrahedron_10; Mesh mesh(spatial_dimension); mesh.read("interpolation.msh"); const ElementType type_facet = mesh.getFacetType(type); Mesh & mesh_facets = mesh.initMeshFacets("mesh_facets"); MeshUtils::buildAllFacets(mesh, mesh_facets); SolidMechanicsModel model(mesh); /// model initialization model.initFull(); Array<Real> & position = mesh.getNodes(); UInt nb_facet = mesh_facets.getNbElement(type_facet); UInt nb_element = mesh.getNbElement(type); /// compute quadrature points positions on facets typedef FEEngineTemplate<IntegratorGauss, ShapeLagrange> MyFEEngineType; model.registerFEEngineObject<MyFEEngineType>("FacetsFEEngine", mesh_facets, spatial_dimension - 1); model.getFEEngine("FacetsFEEngine").initShapeFunctions(); UInt nb_quad_per_facet = model.getFEEngine("FacetsFEEngine").getNbIntegrationPoints(type_facet); UInt nb_tot_quad = nb_quad_per_facet * nb_facet; Array<Real> quad_facets(nb_tot_quad, spatial_dimension); model.getFEEngine("FacetsFEEngine") .interpolateOnIntegrationPoints(position, quad_facets, spatial_dimension, type_facet); Array<Element> & facet_to_element = mesh_facets.getSubelementToElement(type); UInt nb_facet_per_elem = facet_to_element.getNbComponent(); ElementTypeMapArray<Real> element_quad_facet; element_quad_facet.alloc(nb_element * nb_facet_per_elem * nb_quad_per_facet, spatial_dimension, type); ElementTypeMapArray<Real> interpolated_stress("interpolated_stress", ""); interpolated_stress.initialize( mesh, _nb_component = spatial_dimension * spatial_dimension, _spatial_dimension = spatial_dimension); Array<Real> & interp_stress = interpolated_stress(type); interp_stress.resize(nb_element * nb_facet_per_elem * nb_quad_per_facet); Array<Real> & el_q_facet = element_quad_facet(type); for (UInt el = 0; el < nb_element; ++el) { for (UInt f = 0; f < nb_facet_per_elem; ++f) { UInt global_facet = facet_to_element(el, f).element; 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_facets(global_facet * nb_quad_per_facet + q, s); } } } } /// compute quadrature points position of the elements UInt nb_quad_per_element = model.getFEEngine().getNbIntegrationPoints(type); UInt nb_tot_quad_el = nb_quad_per_element * nb_element; Array<Real> quad_elements(nb_tot_quad_el, spatial_dimension); model.getFEEngine().interpolateOnIntegrationPoints(position, quad_elements, spatial_dimension, type); /// assign some values to stresses Array<Real> & stress = const_cast<Array<Real> &>(model.getMaterial(0).getStress(type)); for (UInt q = 0; q < nb_tot_quad_el; ++q) { for (UInt s = 0; s < spatial_dimension * spatial_dimension; ++s) { stress(q, s) = s * function(quad_elements(q, 0), quad_elements(q, 1), quad_elements(q, 2)); } } /// interpolate stresses on facets' quadrature points model.getMaterial(0).initElementalFieldInterpolation(element_quad_facet); model.getMaterial(0).interpolateStress(interpolated_stress); Real tolerance = 1.e-10; /// check results for (UInt el = 0; el < nb_element; ++el) { for (UInt f = 0; f < nb_facet_per_elem; ++f) { for (UInt q = 0; q < nb_quad_per_facet; ++q) { for (UInt s = 0; s < spatial_dimension * spatial_dimension; ++s) { Real x = el_q_facet(el * nb_facet_per_elem * nb_quad_per_facet + f * nb_quad_per_facet + q, 0); Real y = el_q_facet(el * nb_facet_per_elem * nb_quad_per_facet + f * nb_quad_per_facet + q, 1); Real z = el_q_facet(el * nb_facet_per_elem * nb_quad_per_facet + f * nb_quad_per_facet + q, 2); Real theoretical = s * function(x, y, z); Real numerical = interp_stress(el * nb_facet_per_elem * nb_quad_per_facet + f * nb_quad_per_facet + q, s); if (std::abs(theoretical - numerical) > tolerance) { std::cout << "Theoretical and numerical values aren't coincident!" << std::endl; return EXIT_FAILURE; } } } } } std::cout << "OK: Stress interpolation test passed." << std::endl; return EXIT_SUCCESS; }