diff --git a/doc/manual/manual-gettingstarted.tex b/doc/manual/manual-gettingstarted.tex index e450fe238..f56450a04 100644 --- a/doc/manual/manual-gettingstarted.tex +++ b/doc/manual/manual-gettingstarted.tex @@ -1,198 +1,198 @@ \chapter{Getting started} \section{Downloading the code} The \akantu source code can be requested using the form accessible at the URL \url{http://lsms.epfl.ch/akantu}. There, you will be asked to accept the LGPL license terms. \section{Compiling Akantu} Akantu is a \code{cmake} project, so to configure it you can either follow the usual way: \begin{command} > cd akantu > mkdir build > cd build > ccmake [ Set the options that you need ] > make > make install \end{command} Or, you can use a tool we added to help you. You can just use the given \code{Makefile} that handles the \code{cmake} configuration \begin{command} > cd akantu > make config > make > make install \end{command} \section{Writing a \code{main} function\label{sect:common:main}} First of all, \akantu needs to be initialized. There is a memory management included in the core library which allows a correct allocation and de-allocation of vectors, structures and/or objects. Moreover, in parallel computations the initialization procedure can perform the communication setup. This is achieved by a pair of functions (\code{initialize} and \code{finalize}) that are used as follows: \begin{cpp} int main(int argc, char *argv[]) { akantu::initialize(argc, argv); // your code ... akantu::finalize(); } \end{cpp} The \code{initialize} function takes the program parameters which can be interpreted by \akantu in due form. \section{Creating and loading a mesh\label{sect:common:mesh}} \akantu supports meshes generated with Gmsh~\cite{gmsh}, a free software available at \url{http://geuz.org/gmsh/} where a detailed documentation can be found. Consequently, this manual will not provide Gmsh usage directions. Gmsh outputs meshes in \code{.msh} format that can be read by \akantu. In order to import a mesh, it is necessary to create a \code{Mesh} object through the following function calls: \begin{cpp} UInt spatial_dimension = 2; Mesh mesh(spatial_dimension); \end{cpp} The only parameter that has to be specified by the user is the spatial dimension of the problem. Now it is possible to read a \code{.msh} file with a \code{MeshIOMSH} object that takes care of loading a mesh to memory. This step is carried out by: \begin{cpp} MeshIOMSH mesh_io; mesh_io.read("square.msh", mesh); \end{cpp} where the \code{MeshIOMSH} object is first created before being used to read the \code{.msh} file. The mesh file name as well as the \code{Mesh} object must be specified by the user. The \code{MeshIOMSH} object can also write mesh files. This feature is useful to save a mesh that has been modified during a simulation. The \code{write} method takes care of it: \begin{cpp} mesh_io.write("square_modified.msh", mesh); \end{cpp} which works exactly like the \code{read} method. \akantu supports also meshes generated by DIANA (\url{http://tnodiana.com}), but only in reading mode. A similar procedure applies where the only difference is that the \code{MeshIODiana} object should be used instead of the \code{MeshIOMSH} one. Additional mesh readers can be introduced into \akantu by coding new \code{MeshIO} classes. \section{Using \code{Vectors}} Data in \akantu can be stored in data containers implemented by the \code{Vector} class. In its most basic usage the \code{Vector} class implemented in \akantu is similar to the \code{vector} class of the Standard Template Library (STL) for C++. A simple \code{Vector} containing a sequence of \code{nb\_element} values can be generated with: \begin{cpp} Vector example_vector(nb_element); \end{cpp} where \code{type} usually is \code{Real}, \code{UInt} or \code{bool}. Each value is associated to an index, so that data can be accessed by typing: \begin{cpp} type & val = example_vector(index) \end{cpp} \code{Vectors} can also contain a sequence of values for each index. In that case the number of components of each sequence must be specified at the \code{Vector} creation. For example if we want to create a \code{Vector} to store the coordinates (sequences of three values) of ten nodes, the appropriate code is the following: \begin{cpp} UInt nb_nodes = 10; UInt spatial_dimension = 3; Vector position(nb_nodes, spatial_dimension); \end{cpp} In this case the $x$ position of node number 8 will be given by \code{position(7, 0)} (in C++, numbering starts at 0 and not 1). If the number of components for the sequences is not specified, the default value of 1 is used. It is very common in \akantu to loop over vectors to perform a specific treatment. This ranges from geometric calculation on nodal quantities to tensor algebra (in constitutive laws for example). The \code{Vector} object has the possibility to request iterators in order to make the writing of loops easier and enhance readability. For instance, a loop over the nodal coordinates can be performed like: \begin{cpp} //accessing the nodal coordinates Vector Vector nodes = mesh.getNodes(); //creating the iterators Vector::iterator it = nodes.begin(spatial_dimension); Vector::iterator end = nodes.end(spatial_dimension); for (; it != end; ++it){ RVector & coords = (*it); //do what you need .... } \end{cpp} In that example, each \code{RVector} is a geometrical vector of size \code{spatial\_dimension} and the iteration is conveniently performed by the \code{Vector} iterator. The \code{Vector} object is intensively used to store tensor values. In that case it should be specified that the returned object type is a matrix when constructing the iterator. This is done when calling the \code{begin} function. For instance, assuming that we have a \code{Vector} storing stresses, we can loop over the stored tensors by: \begin{cpp} //creating the iterators Vector::iterator it = stresses.begin(spatial_dimension,spatial_dimension); Vector::iterator end = stresses.end(spatial_dimension,spatial_dimension); for (; it != end; ++it){ Matrix & stress = (*it); //do what you need .... } \end{cpp} In that last example the Matrix objects are \code{spatial\_dimension} $\times$ \code{spatial\_dimension} matrices. The light objects \code{Matrix} and \code{RVector} can be used and combined to do most common linear algebra. In general, a mesh consists of several kinds of elements. Consequently the amount of data to be stored can differ for each element type. The straightforward example is the connectivity array, namely the sequences of nodes belonging to each element. In order to easily manage this kind of data, a particular data structure called \code{ByElementTypeVector} is available. It is just a group of \code{Vectors}, each associated to an element type. The following code can retrieve the \code{ByElementTypeVector} which stores the connectivity arrays for a mesh: \begin{cpp} ByElementTypeVector & connectivities = mesh.getConnectivities(); \end{cpp} Then the specific vector associated to a given element type can be obtained by \begin{cpp} Vector & connectivity_triangle = connectivities[_triangle_3]; \end{cpp} where the first order 3-node triangular element was used in the presented piece of code. -%%% Local Variables: +%%% Local Variables: %%% mode: latex %%% TeX-master: "manual" -%%% End: +%%% End: diff --git a/src/fem/element_classes/element_class_quadrangle_4_inline_impl.cc b/src/fem/element_classes/element_class_quadrangle_4_inline_impl.cc index da22a9343..47754b13d 100644 --- a/src/fem/element_classes/element_class_quadrangle_4_inline_impl.cc +++ b/src/fem/element_classes/element_class_quadrangle_4_inline_impl.cc @@ -1,147 +1,158 @@ /** * @file element_class_quadrangle_4_inline_impl.cc * * @author Nicolas Richart * * @date Mon Dec 13 10:45:15 2010 * * @brief Specialization of the element_class class for the type _quadrangle_4 * * @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 . * * @section DESCRIPTION * * @verbatim \eta ^ (-1,1) | (1,1) x---------x | | | | | | --|---------|-----> \xi | | | | | | x---------x (-1,-1) | (1,-1) @endverbatim * * @subsection shapes Shape functions * @f[ * \begin{array}{lll} * N1 = (1 - \xi) (1 - \eta) / 4 * & \frac{\partial N1}{\partial \xi} = - (1 - \eta) / 4 * & \frac{\partial N1}{\partial \eta} = - (1 - \xi) / 4 \\ * N2 = (1 + \xi) (1 - \eta) / 4 \\ * & \frac{\partial N2}{\partial \xi} = (1 - \eta) / 4 * & \frac{\partial N2}{\partial \eta} = - (1 + \xi) / 4 \\ * N3 = (1 + \xi) (1 + \eta) / 4 \\ * & \frac{\partial N3}{\partial \xi} = (1 + \eta) / 4 * & \frac{\partial N3}{\partial \eta} = (1 + \xi) / 4 \\ * N4 = (1 - \xi) (1 + \eta) / 4 * & \frac{\partial N4}{\partial \xi} = - (1 + \eta) / 4 * & \frac{\partial N4}{\partial \eta} = (1 - \xi) / 4 \\ * \end{array} * @f] * * @subsection quad_points Position of quadrature points * @f{eqnarray*}{ * \xi_{q0} &=& 0 \qquad \eta_{q0} = 0 * @f} */ /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_quadrangle_4>::nb_nodes_per_element; template<> UInt ElementClass<_quadrangle_4>::nb_quadrature_points; template<> UInt ElementClass<_quadrangle_4>::spatial_dimension; /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_quadrangle_4>::computeShapes(const Real * natural_coords, Real * shapes) { /// Natural coordinates const Real * c = natural_coords; shapes[0] = .25 * (1 - c[0]) * (1 - c[1]); /// N1(q_0) shapes[1] = .25 * (1 + c[0]) * (1 - c[1]); /// N2(q_0) shapes[2] = .25 * (1 + c[0]) * (1 + c[1]); /// N3(q_0) shapes[3] = .25 * (1 - c[0]) * (1 + c[1]); /// N4(q_0) } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_quadrangle_4>::computeDNDS(const Real * natural_coords, Real * dnds) { /** * @f[ * dnds = \left( * \begin{array}{cccc} * \frac{\partial N1}{\partial \xi} & \frac{\partial N2}{\partial \xi} * & \frac{\partial N3}{\partial \xi} & \frac{\partial N4}{\partial \xi}\\ * \frac{\partial N1}{\partial \eta} & \frac{\partial N2}{\partial \eta} * & \frac{\partial N3}{\partial \eta} & \frac{\partial N4}{\partial \eta} * \end{array} * \right) * @f] */ const Real * c = natural_coords; dnds[0] = - .25 * (1 - c[1]); dnds[1] = .25 * (1 - c[1]); dnds[2] = .25 * (1 + c[1]); dnds[3] = - .25 * (1 + c[1]); dnds[4] = - .25 * (1 - c[0]); dnds[5] = - .25 * (1 + c[0]); dnds[6] = .25 * (1 + c[0]); dnds[7] = .25 * (1 - c[0]); } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_quadrangle_4>::computeJacobian(const Real * dxds, const UInt dimension, Real & jac){ if (dimension == spatial_dimension){ Real det_dxds = Math::det2(dxds); jac = det_dxds; } else { AKANTU_DEBUG_TO_IMPLEMENT(); } } /* -------------------------------------------------------------------------- */ template<> inline Real ElementClass<_quadrangle_4>::getInradius(const Real * coord) { Real a = Math::distance_2d(coord + 0, coord + 2); Real b = Math::distance_2d(coord + 2, coord + 4); Real c = Math::distance_2d(coord + 4, coord + 6); Real d = Math::distance_2d(coord + 6, coord + 0); // Real septimetre = (a + b + c + d) / 2.; // Real p = Math::distance_2d(coord + 0, coord + 4); // Real q = Math::distance_2d(coord + 2, coord + 6); // Real area = sqrt(4*(p*p * q*q) - (a*a + b*b + c*c + d*d)*(a*a + c*c - b*b - d*d)) / 4.; // Real h = sqrt(area); // to get a length // Real h = area / septimetre; // formula of inradius for circumscritable quadrelateral Real h = std::min(a, std::min(b, std::min(c, d))); return h; } + +/* -------------------------------------------------------------------------- */ + +template<> inline bool ElementClass<_quadrangle_4>::contains(const types::RVector & natural_coords) { + if (natural_coords[0] < -1.) return false; + if (natural_coords[0] > 1.) return false; + if (natural_coords[1] < -1.) return false; + if (natural_coords[1] > 1.) return false; + return true; +} +/* -------------------------------------------------------------------------- */ diff --git a/src/fem/fem_template.hh b/src/fem/fem_template.hh index 618be0e70..efaa09f02 100644 --- a/src/fem/fem_template.hh +++ b/src/fem/fem_template.hh @@ -1,272 +1,272 @@ /** * @file fem_template.hh * * @author Guillaume Anciaux * * @date Tue Feb 15 16:32:44 2011 * * @brief templated class that calls integration and shape objects * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ #ifndef __AKANTU_FEM_TEMPLATE_HH__ #define __AKANTU_FEM_TEMPLATE_HH__ /* -------------------------------------------------------------------------- */ #include "fem.hh" #include "integrator.hh" #include "shape_functions.hh" #include "shape_lagrange.hh" #include "shape_linked.hh" #include "integrator_gauss.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template class FEMTemplate : public FEM{ /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: FEMTemplate(Mesh & mesh, UInt spatial_dimension = 0, ID id = "fem", MemoryID memory_id = 0); virtual ~FEMTemplate(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// pre-compute all the shape functions, their derivatives and the jacobians void initShapeFunctions(const GhostType & ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ /* Integration method bridges */ /* ------------------------------------------------------------------------ */ - /// integrate f for all elements of type "type" void integrate(const Vector & f, Vector &intf, UInt nb_degree_of_freedom, const ElementType & type, const GhostType & ghost_type = _not_ghost, const Vector * filter_elements = NULL) const; /// integrate a scalar value on all elements of type "type" Real integrate(const Vector & f, const ElementType & type, const GhostType & ghost_type = _not_ghost, const Vector * filter_elements = NULL) const; /// integrate one element scalar value on all elements of type "type" virtual Real integrate(const types::RVector & f, const ElementType & type, UInt index, const GhostType & ghost_type = _not_ghost) const; /// integrate partially around a quadrature point (@f$ intf_q = f_q * J_q * w_q @f$) void integrateOnQuadraturePoints(const Vector & f, Vector &intf, UInt nb_degree_of_freedom, const ElementType & type, const GhostType & ghost_type = _not_ghost, const Vector * filter_elements = NULL) const; /// get the number of quadrature points UInt getNbQuadraturePoints(const ElementType & type, const GhostType & ghost_type = _not_ghost) const; /// get shapes precomputed const Vector & getShapes(const ElementType & type, const GhostType & ghost_type = _not_ghost) const; /// get the derivatives of shapes const Vector & getShapesDerivatives(const ElementType & type, const GhostType & ghost_type = _not_ghost, UInt id=0) const; /// get quadrature points const Vector & getQuadraturePoints(const ElementType & type, const GhostType & ghost_type = _not_ghost) const; /* ------------------------------------------------------------------------ */ /* Shape method bridges */ /* ------------------------------------------------------------------------ */ - + /// compute the gradient of a nodal field on the quadrature points void gradientOnQuadraturePoints(const Vector &u, Vector &nablauq, const UInt nb_degree_of_freedom, const ElementType & type, const GhostType & ghost_type = _not_ghost, const Vector * filter_elements = NULL) const; + /// interpolate a nodal field on the quadrature points void interpolateOnQuadraturePoints(const Vector &u, Vector &uq, UInt nb_degree_of_freedom, const ElementType & type, const GhostType & ghost_type = _not_ghost, const Vector * filter_elements = NULL) const; /// find natural coords from real coords provided an element void inverseMap(const types::RVector & real_coords, UInt element, const ElementType & type, types::RVector & 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 types::RVector & 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 types::RVector & real_coords, UInt element, const ElementType & type, types::RVector & shapes, const GhostType & ghost_type = _not_ghost) const; /* ------------------------------------------------------------------------ */ /* Other methods */ /* ------------------------------------------------------------------------ */ /// pre-compute normals on control points void computeNormalsOnControlPoints(const GhostType & ghost_type = _not_ghost); void computeNormalsOnControlPoints(const Vector & field, const GhostType & ghost_type = _not_ghost); void computeNormalsOnControlPoints(const Vector & field, Vector & normal, const ElementType & type, const GhostType & ghost_type = _not_ghost) const; /// function to print the contain of the class // virtual void printself(std::ostream & stream, int indent = 0) const{}; void assembleFieldLumped(const Vector & field_1, UInt nb_degree_of_freedom, Vector & lumped, const Vector & equation_number, ElementType type, const GhostType & ghost_type = _not_ghost) const; void assembleFieldMatrix(const Vector & field, UInt nb_degree_of_freedom, SparseMatrix & matrix, ElementType type, const GhostType & ghost_type = _not_ghost) const; private: template void assembleLumpedTemplate(const Vector & field_1, UInt nb_degree_of_freedom, Vector & lumped, const Vector & equation_number, const GhostType & ghost_type) const; /// @f$ \tilde{M}_{i} = \sum_j M_{ij} = \sum_j \int \rho \varphi_i \varphi_j dV = \int \rho \varphi_i dV @f$ template void assembleLumpedRowSum(const Vector & field_1, UInt nb_degree_of_freedom, Vector & lumped, const Vector & equation_number, const GhostType & ghost_type) const; /// @f$ \tilde{M}_{i} = c * M_{ii} = \int_{V_e} \rho dV @f$ template void assembleLumpedDiagonalScaling(const Vector & field_1, UInt nb_degree_of_freedom, Vector & lumped, const Vector & equation_number, const GhostType & ghost_type) const; template void assembleFieldMatrix(const Vector & field, UInt nb_degree_of_freedom, SparseMatrix & matrix, const GhostType & ghost_type) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: const ShapeFunctions & getShapeFunctionsInterface() const { return shape_functions; }; const Shape & getShapeFunctions() const { return shape_functions; }; const Integrator & getIntegratorInterface() const { return integrator; }; const Integ & getIntegrator() const { return integrator; }; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: Integ integrator; Shape shape_functions; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "fem_template_inline_impl.cc" #include "fem_template_tmpl.hh" /* -------------------------------------------------------------------------- */ #if defined(AKANTU_COHESIVE_ELEMENT) # include "fem_template_cohesive_inline_impl.cc" # include "fem_template_cohesive_tmpl.hh" #endif /// standard output stream operator // inline std::ostream & operator <<(std::ostream & stream, const FEMTemplate & _this) // { // _this.printself(stream); // return stream; // } __END_AKANTU__ // #if defined(AKANTU_COHESIVE_ELEMENT) // # include "shape_cohesive.hh" // # include "integrator_cohesive.hh" // #endif #endif /* __AKANTU_FEM_TEMPLATE_HH__ */ diff --git a/src/fem/mesh.hh b/src/fem/mesh.hh index 159294fc9..658278e87 100644 --- a/src/fem/mesh.hh +++ b/src/fem/mesh.hh @@ -1,603 +1,603 @@ /** * @file mesh.hh * * @author Guillaume Anciaux * @author Marco Vocialta * @author Nicolas Richart * * @date Fri Jun 18 11:47:19 2010 * * @brief the class representing the meshes * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_HH__ #define __AKANTU_MESH_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "aka_vector.hh" #include "element_class.hh" #include "by_element_type.hh" #include "aka_event_handler.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* Element */ /* -------------------------------------------------------------------------- */ class Element { public: Element(ElementType type = _not_defined, UInt element = 0, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) : type(type), element(element), ghost_type(ghost_type), kind(kind) {}; Element(const Element & element) { this->type = element.type; this->element = element.element; this->ghost_type = element.ghost_type; this->kind = element.kind; } inline bool operator==(const Element & elem) const { return ((element == elem.element) && (type == elem.type) && (ghost_type == elem.ghost_type) && (kind == elem.kind)); } inline bool operator!=(const Element & elem) const { return ((element != elem.element) || (type != elem.type) || (ghost_type != elem.ghost_type) || (kind != elem.kind)); } bool operator<(const Element& rhs) const { bool res = ((this->kind < rhs.kind) || ((this->kind == rhs.kind) && ((this->ghost_type < rhs.ghost_type) || ((this->ghost_type == rhs.ghost_type) && ((this->type < rhs.type) || ((this->type == rhs.type) && (this->element < rhs.element))))))); return res; } virtual ~Element() {}; /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; public: ElementType type; UInt element; GhostType ghost_type; ElementKind kind; }; struct CompElementLess { bool operator() (const Element& lhs, const Element& rhs) const { return lhs < rhs; } }; extern const Element ElementNull; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* Mesh modifications events */ /* -------------------------------------------------------------------------- */ template class MeshEvent { public: const Vector & getList() const { return list; } Vector & getList() { return list; } protected: Vector list; }; class Mesh; class NewNodesEvent : public MeshEvent { }; class RemovedNodesEvent : public MeshEvent { public: inline RemovedNodesEvent(const Mesh & mesh); AKANTU_GET_MACRO_NOT_CONST(NewNumbering, new_numbering, Vector &); AKANTU_GET_MACRO(NewNumbering, new_numbering, const Vector &); private: Vector new_numbering; }; class NewElementsEvent : public MeshEvent { }; class RemovedElementsEvent : public MeshEvent { public: inline RemovedElementsEvent(const Mesh & mesh); AKANTU_GET_MACRO(NewNumbering, new_numbering, const ByElementTypeUInt &); AKANTU_GET_MACRO_NOT_CONST(NewNumbering, new_numbering, ByElementTypeUInt &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(NewNumbering, new_numbering, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(NewNumbering, new_numbering, UInt); protected: ByElementTypeUInt new_numbering; }; /* -------------------------------------------------------------------------- */ class MeshEventHandler { public: virtual ~MeshEventHandler() {}; /* ------------------------------------------------------------------------ */ /* Internal code */ /* ------------------------------------------------------------------------ */ protected: inline void sendEvent(const NewNodesEvent & event) { onNodesAdded (event.getList()); } inline void sendEvent(const RemovedNodesEvent & event) { onNodesRemoved(event.getList(), event.getNewNumbering()); } inline void sendEvent(const NewElementsEvent & event) { onElementsAdded (event.getList()); } inline void sendEvent(const RemovedElementsEvent & event) { onElementsRemoved(event.getList(), event.getNewNumbering()); } template friend class EventHandlerManager; /* ------------------------------------------------------------------------ */ /* Interface */ /* ------------------------------------------------------------------------ */ public: virtual void onNodesAdded (__attribute__((unused)) const Vector & nodes_list) { } virtual void onNodesRemoved(__attribute__((unused)) const Vector & nodes_list, __attribute__((unused)) const Vector & new_numbering) { } virtual void onElementsAdded (__attribute__((unused)) const Vector & elements_list) { } virtual void onElementsRemoved(__attribute__((unused)) const Vector & elements_list, __attribute__((unused)) const ByElementTypeUInt & new_numbering) { } }; /* -------------------------------------------------------------------------- */ /* Mesh */ /* -------------------------------------------------------------------------- */ /** * @class Mesh this contain the coordinates of the nodes in the Mesh.nodes * Vector, and the connectivity. The connectivity are stored in by element * types. * * To know all the element types present in a mesh you can get the * Mesh::ConnectivityTypeList * * In order to loop on all element you have to loop on all types like this : * @code Mesh::type_iterator it = mesh.firstType(dim, ghost_type); Mesh::type_iterator end = mesh.lastType(dim, ghost_type); for(; it != end; ++it) { UInt nb_element = mesh.getNbElement(*it); const Vector & conn = mesh.getConnectivity(*it); for(UInt e = 0; e < nb_element; ++e) { ... } } @endcode */ class Mesh : protected Memory, public EventHandlerManager { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /// constructor that create nodes coordinates array Mesh(UInt spatial_dimension, const ID & id = "mesh", const MemoryID & memory_id = 0); /// constructor that use an existing nodes coordinates array, by knowing its ID Mesh(UInt spatial_dimension, const ID & nodes_id, const ID & id = "mesh", const MemoryID & memory_id = 0); /** * constructor that use an existing nodes coordinates * array, by getting the vector of coordinates */ Mesh(UInt spatial_dimension, Vector & nodes, const ID & id = "mesh", const MemoryID & memory_id = 0); virtual ~Mesh(); /// @typedef ConnectivityTypeList list of the types present in a Mesh typedef std::set ConnectivityTypeList; // typedef Vector * NormalsMap[_max_element_type]; private: /// initialize the connectivity to NULL and other stuff void init(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; /// function that computes the bounding box (fills xmin, xmax) void computeBoundingBox(); /// init a by-element-type real vector with provided ids template void initByElementTypeVector(ByElementTypeVector & v, UInt nb_component, UInt spatial_dimension, const bool & flag_nb_node_per_elem_multiply = false, ElementKind element_kind = _ek_regular, bool size_to_nb_element = false) const; /// @todo: think about nicer way to do it /// extract coordinates of nodes from an element template inline void extractNodalValuesFromElement(const Vector & nodal_values, T * elemental_values, UInt * connectivity, UInt n_nodes, UInt nb_degree_of_freedom) const; /// extract coordinates of nodes from a reversed element inline void extractNodalCoordinatesFromPBCElement(Real * local_coords, UInt * connectivity, UInt n_nodes); /// convert a element to a linearized element - inline UInt elementToLinearized(const Element & elem); + inline UInt elementToLinearized(const Element & elem) const; /// convert a linearized element to an element - inline Element linearizedToElement (UInt linearized_element); + inline Element linearizedToElement (UInt linearized_element) const; /// update the types offsets array for the conversions inline void updateTypesOffsets(const GhostType & ghost_type); /// add a Vector of connectivity for the type . inline void addConnectivityType(const ElementType & type); /* ------------------------------------------------------------------------ */ template inline void sendEvent(Event & event) { if(event.getList().getSize() != 0) EventHandlerManager::sendEvent(event); } /* ------------------------------------------------------------------------ */ template inline void removeNodesFromVector(Vector & vect, const Vector & new_numbering); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(ID, id, const ID &); /// get the spatial dimension of the mesh = number of component of the coordinates AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the nodes Vector aka coordinates AKANTU_GET_MACRO(Nodes, *nodes, const Vector &); AKANTU_GET_MACRO_NOT_CONST(Nodes, *nodes, Vector &); /// get the number of nodes AKANTU_GET_MACRO(NbNodes, nodes->getSize(), UInt); /// get the Vector of global ids of the nodes (only used in parallel) AKANTU_GET_MACRO(GlobalNodesIds, *nodes_global_ids, const Vector &); /// get the global id of a node inline UInt getNodeGlobalId(UInt local_id) const; /// get the global number of nodes inline UInt getNbGlobalNodes() const; /// get the nodes type Vector AKANTU_GET_MACRO(NodesType, *nodes_type, const Vector &); inline Int getNodeType(UInt local_id) const; /// say if a node is a pure ghost node inline bool isPureGhostNode(UInt n) const; /// say if a node is pur local or master node inline bool isLocalOrMasterNode(UInt n) const; inline bool isLocalNode(UInt n) const; inline bool isMasterNode(UInt n) const; inline bool isSlaveNode(UInt n) const; AKANTU_GET_MACRO(XMin, lower_bounds[0], Real); AKANTU_GET_MACRO(YMin, lower_bounds[1], Real); AKANTU_GET_MACRO(ZMin, lower_bounds[2], Real); AKANTU_GET_MACRO(XMax, upper_bounds[0], Real); AKANTU_GET_MACRO(YMax, upper_bounds[1], Real); AKANTU_GET_MACRO(ZMax, upper_bounds[2], Real); inline void getLowerBounds(Real * lower) const; inline void getUpperBounds(Real * upper) const; inline void getLocalLowerBounds(Real * lower) const; inline void getLocalUpperBounds(Real * upper) const; /// get the number of surfaces AKANTU_GET_MACRO(NbSurfaces, nb_surfaces, UInt); /// get the connectivity Vector for a given type AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Connectivity, connectivities, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Connectivity, connectivities, UInt); AKANTU_GET_MACRO(Connectivities, connectivities, const ByElementTypeVector &); /// @todo take out this set, if mesh can read surface id /// set the number of surfaces AKANTU_SET_MACRO(NbSurfaces, nb_surfaces, UInt); /// get the number of element of a type in the mesh inline UInt getNbElement(const ElementType & type, const GhostType & ghost_type = _not_ghost) const; // /// get the number of ghost element of a type in the mesh // inline UInt getNbGhostElement(const ElementType & type) const; /// get the connectivity list either for the elements or the ghost elements inline const ConnectivityTypeList & getConnectivityTypeList(const GhostType & ghost_type = _not_ghost) const; /// compute the barycenter of a given element inline void getBarycenter(UInt element, const ElementType & type, Real * barycenter, GhostType ghost_type = _not_ghost) const; /// get the element connected to a subelement AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementToSubelement, element_to_subelement, Vector); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementToSubelement, element_to_subelement, Vector); /// get the subelement connected to an element AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(SubelementToElement, subelement_to_element, Element); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(SubelementToElement, subelement_to_element, Element); /// get the surface values of facets AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(SurfaceID, surface_id, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(SurfaceID, surface_id, UInt); /// set the int data to the surface id vectors void setSurfaceIDsFromIntData(const std::string & data_name); inline const Vector & getUIntData(const ElementType & el_type, const std::string & data_name, const GhostType & ghost_type = _not_ghost) const; /* ------------------------------------------------------------------------ */ /* Wrappers on ElementClass functions */ /* ------------------------------------------------------------------------ */ public: /// get the number of nodes per element for a given element type static inline UInt getNbNodesPerElement(const ElementType & type); /// get the number of nodes per element for a given element type considered as /// a first order element static inline ElementType getP1ElementType(const ElementType & type); /// get the kind of the element type static inline ElementKind getKind(const ElementType & type); /// get spatial dimension of a type of element static inline UInt getSpatialDimension(const ElementType & type); /// get number of facets of a given element type static inline UInt getNbFacetsPerElement(const ElementType & type); /// get local connectivity of a facet for a given facet type static inline UInt ** getFacetLocalConnectivity(const ElementType & type); /// get the type of the surface element associated to a given element static inline ElementType getFacetElementType(const ElementType & type); /// get the pointer to the list of elements for a given type inline Vector * getReversedElementsPBCPointer(const ElementType & type); /* ------------------------------------------------------------------------ */ /* Element type Iterator */ /* ------------------------------------------------------------------------ */ typedef ByElementTypeUInt::type_iterator type_iterator; inline type_iterator firstType(UInt dim = 0, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const { return connectivities.firstType(dim, ghost_type, kind); } inline type_iterator lastType(UInt dim = 0, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const { return connectivities.lastType(dim, ghost_type, kind); } /* ------------------------------------------------------------------------ */ /* Private methods for friends */ /* ------------------------------------------------------------------------ */ private: friend class MeshIOMSH; friend class MeshIOMSHStruct; friend class MeshIODiana; friend class MeshUtils; friend class DistributedSynchronizer; AKANTU_GET_MACRO(NodesPointer, nodes, Vector *); /// get a pointer to the nodes_global_ids Vector and create it if necessary inline Vector * getNodesGlobalIdsPointer(); /// get a pointer to the nodes_type Vector and create it if necessary inline Vector * getNodesTypePointer(); /// get a pointer to the connectivity Vector for the given type and create it if necessary inline Vector * getConnectivityPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); // inline Vector * getNormalsPointer(ElementType type) const; /// get a pointer to the surface_id Vector for the given type and create it if necessary inline Vector * getSurfaceIDPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); /// get the UIntDataMap for a given ElementType inline UIntDataMap & getUIntDataMap(const ElementType & el_type, const GhostType & ghost_type = _not_ghost); /// get the IntDataMap pointer (modifyable) for a given ElementType inline Vector * getUIntDataPointer(const ElementType & el_type, const std::string & data_name, const GhostType & ghost_type = _not_ghost); /// get a pointer to the element_to_subelement Vector for the given type and create it if necessary inline Vector > * getElementToSubelementPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); /// get a pointer to the subelement_to_element Vector for the given type and create it if necessary inline Vector * getSubelementToElementPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// id of the mesh ID id; /// array of the nodes coordinates Vector * nodes; /// global node ids Vector * nodes_global_ids; /// node type, -3 pure ghost, -2 master for the node, -1 normal node, i in /// [0-N] slave node and master is proc i Vector * nodes_type; /// global number of nodes; UInt nb_global_nodes; /// boolean to know if the nodes have to be deleted with the mesh or not bool created_nodes; /// all class of elements present in this mesh (for heterogenous meshes) ByElementTypeUInt connectivities; /// map to normals for all class of elements present in this mesh ByElementTypeReal normals; /// list of all existing types in the mesh ConnectivityTypeList type_set; /// the spatial dimension of this mesh UInt spatial_dimension; /// types offsets Vector types_offsets; /// list of all existing types in the mesh ConnectivityTypeList ghost_type_set; /// ghost types offsets Vector ghost_types_offsets; /// number of surfaces present in this mesh UInt nb_surfaces; /// surface id of the surface elements in this mesh ByElementTypeUInt surface_id; /// min of coordinates Real lower_bounds[3]; /// max of coordinates Real upper_bounds[3]; /// size covered by the mesh on each direction Real size[3]; /// local min of coordinates Real local_lower_bounds[3]; /// local max of coordinates Real local_upper_bounds[3]; /// List of elements connected to subelements ByElementTypeVector > element_to_subelement; /// List of subelements connected to elements ByElementTypeVector subelement_to_element; // /// list of elements that are reversed due to pbc // ByElementTypeUInt reversed_elements_pbc; // /// direction in which pbc are to be applied // UInt pbc_directions[3]; /// list of the vectors corresponding to tags in the mesh ByElementTypeUIntDataMap uint_data; }; /* -------------------------------------------------------------------------- */ /* Inline functions */ /* -------------------------------------------------------------------------- */ /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Element & _this) { _this.printself(stream); return stream; } #if defined (AKANTU_INCLUDE_INLINE_IMPL) # include "mesh_inline_impl.cc" #endif #include "by_element_type_tmpl.hh" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Mesh & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_MESH_HH__ */ diff --git a/src/fem/mesh_inline_impl.cc b/src/fem/mesh_inline_impl.cc index bfa1ad0e8..8e03c0f2b 100644 --- a/src/fem/mesh_inline_impl.cc +++ b/src/fem/mesh_inline_impl.cc @@ -1,553 +1,550 @@ /** * @file mesh_inline_impl.cc * * @author Guillaume Anciaux * @author Marco Vocialta * @author Nicolas Richart * * @date Thu Jul 15 00:41:12 2010 * * @brief Implementation of the inline functions of the mesh class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ __END_AKANTU__ #if defined(AKANTU_COHESIVE_ELEMENT) # include "cohesive_element.hh" #endif __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ inline RemovedNodesEvent::RemovedNodesEvent(const Mesh & mesh) : new_numbering(mesh.getNbNodes(), 1, "new_numbering") { } /* -------------------------------------------------------------------------- */ inline RemovedElementsEvent::RemovedElementsEvent(const Mesh & mesh) : new_numbering("new_numbering", mesh.getID()) { } /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent(RemovedElementsEvent & event) { - if(event.getList().getSize() != 0) { - connectivities.onElementsRemoved(event.getNewNumbering()); - - EventHandlerManager::sendEvent(event); - } + connectivities.onElementsRemoved(event.getNewNumbering()); + EventHandlerManager::sendEvent(event); } /* -------------------------------------------------------------------------- */ template <> inline void Mesh::sendEvent(RemovedNodesEvent & event) { if(created_nodes) removeNodesFromVector(*nodes , event.getNewNumbering()); if(nodes_global_ids) removeNodesFromVector(*nodes_global_ids, event.getNewNumbering()); if(nodes_type) removeNodesFromVector(*nodes_type , event.getNewNumbering()); EventHandlerManager::sendEvent(event); } /* -------------------------------------------------------------------------- */ template inline void Mesh::removeNodesFromVector(Vector & vect, const Vector & new_numbering) { Vector 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)) { memcpy(tmp.storage() + new_i * nb_component, vect.storage() + i * nb_component, nb_component * sizeof(T)); ++new_nb_nodes; } } tmp.resize(new_nb_nodes); vect.copy(tmp); } /* -------------------------------------------------------------------------- */ -inline UInt Mesh::elementToLinearized(const Element & elem) { +inline UInt Mesh::elementToLinearized(const Element & elem) const { AKANTU_DEBUG_ASSERT(elem.type < _max_element_type && elem.element < types_offsets.values[elem.type+1], "The element " << elem << "does not exists in the mesh " << id); return types_offsets.values[elem.type] + elem.element; } /* -------------------------------------------------------------------------- */ -inline Element Mesh::linearizedToElement (UInt linearized_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(t != _max_element_type, "The linearized element " << linearized_element << "does not exists in the mesh " << id); --t; return Element(ElementType(t), linearized_element - types_offsets.values[t]); } /* -------------------------------------------------------------------------- */ inline void Mesh::updateTypesOffsets(const GhostType & ghost_type) { types_offsets.clear(); type_iterator it = firstType(0, ghost_type); type_iterator last = lastType(0, ghost_type); for (; it != last; ++it) types_offsets(*it) = connectivities(*it, 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 Vector * Mesh::getNodesGlobalIdsPointer() { AKANTU_DEBUG_IN(); if(nodes_global_ids == NULL) { std::stringstream sstr; sstr << id << ":nodes_global_ids"; nodes_global_ids = &(alloc(sstr.str(), nodes->getSize(), 1)); } AKANTU_DEBUG_OUT(); return nodes_global_ids; } /* -------------------------------------------------------------------------- */ inline Vector * Mesh::getNodesTypePointer() { AKANTU_DEBUG_IN(); if(nodes_type == NULL) { std::stringstream sstr; sstr << id << ":nodes_type"; nodes_type = &(alloc(sstr.str(), nodes->getSize(), 1, -1)); } AKANTU_DEBUG_OUT(); return nodes_type; } /* -------------------------------------------------------------------------- */ inline Vector * Mesh::getConnectivityPointer(const ElementType & type, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); Vector * 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 Vector * Mesh::getSurfaceIDPointer(const ElementType & type, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); Vector * tmp; if(!surface_id.exists(type, ghost_type)) { tmp = &(surface_id.alloc(0, 1, type, ghost_type)); AKANTU_DEBUG_INFO("The surface id vector for the type " << type << " created"); } else { tmp = &(surface_id(type, ghost_type)); } AKANTU_DEBUG_OUT(); return tmp; } /* -------------------------------------------------------------------------- */ inline Vector > * Mesh::getElementToSubelementPointer(const ElementType & type, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); Vector > * tmp; if(!element_to_subelement.exists(type, ghost_type)) { tmp = &(element_to_subelement.alloc(0, 1, type, ghost_type)); AKANTU_DEBUG_INFO("The element_to_subelement vector for the type " << type << " created"); } else { tmp = &(element_to_subelement(type, ghost_type)); } AKANTU_DEBUG_OUT(); return tmp; } /* -------------------------------------------------------------------------- */ inline Vector * Mesh::getSubelementToElementPointer(const ElementType & type, const GhostType & ghost_type) { AKANTU_DEBUG_IN(); Vector * tmp; if(!subelement_to_element.exists(type, ghost_type)) { UInt nb_facets = getNbFacetsPerElement(type); tmp = &(subelement_to_element.alloc(0, nb_facets, type, ghost_type)); AKANTU_DEBUG_INFO("The subelement_to_element vector for the type " << type << " created"); } else { tmp = &(subelement_to_element(type, ghost_type)); } AKANTU_DEBUG_OUT(); return tmp; } /* -------------------------------------------------------------------------- */ inline Vector * Mesh::getUIntDataPointer(const ElementType & el_type, const std::string & data_name, const GhostType & ghost_type) { // AKANTU_DEBUG_IN(); Vector * data; // if(!uint_data.exists(el_type, ghost_type)){ // uint_data(UIntDataMap(), el_type, ghost_type); // } UIntDataMap & map = uint_data(el_type, ghost_type); UIntDataMap::iterator it = map.find(data_name); if(it == map.end()) { data = new Vector(0, 1, data_name); map[data_name] = data; } else { data = it->second; } // AKANTU_DEBUG_OUT(); return data; } /* -------------------------------------------------------------------------- */ inline const Vector & Mesh::getUIntData(const ElementType & el_type, const std::string & data_name, const GhostType & ghost_type) const { AKANTU_DEBUG_IN(); const UIntDataMap & map = uint_data(el_type, ghost_type); UIntDataMap::const_iterator it = map.find(data_name); AKANTU_DEBUG_ASSERT(it != map.end(), "No data named " << data_name << " in the mesh " << id); AKANTU_DEBUG_OUT(); return *(it->second); } /* -------------------------------------------------------------------------- */ inline UIntDataMap & Mesh::getUIntDataMap(const ElementType & el_type, const GhostType & ghost_type) { // AKANTU_DEBUG_ASSERT(uint_data.exists(el_type, ghost_type), // "No UIntDataMap for the type (" << ghost_type << ":" << el_type << ")"); return uint_data(el_type, ghost_type); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbElement(const ElementType & type, const GhostType & ghost_type) const { AKANTU_DEBUG_IN(); try { const Vector & conn = connectivities(type, ghost_type); AKANTU_DEBUG_OUT(); return conn.getSize(); } catch (...) { AKANTU_DEBUG_OUT(); return 0; } } /* -------------------------------------------------------------------------- */ inline void Mesh::getBarycenter(UInt element, const ElementType & type, Real * barycenter, GhostType ghost_type) const { AKANTU_DEBUG_IN(); UInt * conn_val = getConnectivity(type, ghost_type).values; UInt nb_nodes_per_element = getNbNodesPerElement(type); Real local_coord[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); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbNodesPerElement(const ElementType & type) { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = 0; #define GET_NB_NODES_PER_ELEMENT(type) \ nb_nodes_per_element = ElementClass::getNbNodesPerElement() #define GET_NB_NODES_PER_ELEMENT_COHESIVE(type) \ nb_nodes_per_element = CohesiveElement::getNbNodesPerElement() AKANTU_BOOST_ELEMENT_SWITCH(GET_NB_NODES_PER_ELEMENT, AKANTU_REGULAR_ELEMENT_TYPE, GET_NB_NODES_PER_ELEMENT_COHESIVE, AKANTU_COHESIVE_ELEMENT_TYPE); #undef GET_NB_NODES_PER_ELEMENT #undef GET_NB_NODES_PER_ELEMENT_COHESIVE AKANTU_DEBUG_OUT(); return nb_nodes_per_element; } /* -------------------------------------------------------------------------- */ inline ElementType Mesh::getP1ElementType(const ElementType & type) { AKANTU_DEBUG_IN(); ElementType element_p1 = _not_defined; #define GET_ELEMENT_P1(type) \ element_p1 = ElementClass::getP1ElementType() AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(GET_ELEMENT_P1); #undef GET_NB_NODES_PER_ELEMENT_P1 AKANTU_DEBUG_OUT(); return element_p1; } /* -------------------------------------------------------------------------- */ inline ElementKind Mesh::getKind(const ElementType & type) { AKANTU_DEBUG_IN(); ElementKind kind = _ek_not_defined; #define GET_KIND(type) \ kind = ElementClass::getKind() AKANTU_BOOST_ELEMENT_SWITCH(GET_KIND, AKANTU_REGULAR_ELEMENT_TYPE, GET_KIND, AKANTU_COHESIVE_ELEMENT_TYPE); #undef GET_KIND AKANTU_DEBUG_OUT(); return kind; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getSpatialDimension(const ElementType & type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = 0; #define GET_SPATIAL_DIMENSION(type) \ spatial_dimension = ElementClass::getSpatialDimension() #define GET_SPATIAL_DIMENSION_COHESIVE(type) \ spatial_dimension = CohesiveElement::getSpatialDimension() AKANTU_BOOST_ELEMENT_SWITCH(GET_SPATIAL_DIMENSION, AKANTU_REGULAR_ELEMENT_TYPE, GET_SPATIAL_DIMENSION_COHESIVE, AKANTU_COHESIVE_ELEMENT_TYPE); #undef GET_SPATIAL_DIMENSION #undef GET_SPATIAL_DIMENSION_COHESIVE AKANTU_DEBUG_OUT(); return spatial_dimension; } /* -------------------------------------------------------------------------- */ inline ElementType Mesh::getFacetElementType(const ElementType & type) { AKANTU_DEBUG_IN(); ElementType surface_type = _not_defined; #define GET_FACET_TYPE(type) \ surface_type = ElementClass::getFacetElementType() #define GET_FACET_TYPE_COHESIVE(type) \ surface_type = CohesiveElement::getFacetElementType() AKANTU_BOOST_ELEMENT_SWITCH(GET_FACET_TYPE, AKANTU_REGULAR_ELEMENT_TYPE, GET_FACET_TYPE_COHESIVE, AKANTU_COHESIVE_ELEMENT_TYPE); #undef GET_FACET_TYPE #undef GET_FACET_TYPE_COHESIVE AKANTU_DEBUG_OUT(); return surface_type; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbFacetsPerElement(const ElementType & type) { AKANTU_DEBUG_IN(); UInt n_facet = 0; #define GET_NB_FACET(type) \ n_facet = ElementClass::getNbFacetsPerElement() #define GET_NB_FACET_COHESIVE(type) \ n_facet = CohesiveElement::getNbFacetsPerElement() AKANTU_BOOST_ELEMENT_SWITCH(GET_NB_FACET, AKANTU_REGULAR_ELEMENT_TYPE, GET_NB_FACET_COHESIVE, AKANTU_COHESIVE_ELEMENT_TYPE); #undef GET_NB_FACET #undef GET_NB_FACET_COHESIVE AKANTU_DEBUG_OUT(); return n_facet; } /* -------------------------------------------------------------------------- */ inline UInt ** Mesh::getFacetLocalConnectivity(const ElementType & type) { AKANTU_DEBUG_IN(); UInt ** facet_conn = NULL; #define GET_FACET_CON(type) \ facet_conn = ElementClass::getFacetLocalConnectivityPerElement() #define GET_FACET_CON_COHESIVE(type) \ facet_conn = CohesiveElement::getFacetLocalConnectivityPerElement() AKANTU_BOOST_ELEMENT_SWITCH(GET_FACET_CON, AKANTU_REGULAR_ELEMENT_TYPE, GET_FACET_CON_COHESIVE, AKANTU_COHESIVE_ELEMENT_TYPE); #undef GET_FACET_CON #undef GET_FACET_CON_COHESIVE AKANTU_DEBUG_OUT(); return facet_conn; } /* -------------------------------------------------------------------------- */ template inline void Mesh::extractNodalValuesFromElement(const Vector & 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)); } } /* -------------------------------------------------------------------------- */ #define DECLARE_GET_BOUND(Var, var) \ inline void Mesh::get##Var##Bounds(Real * var) const { \ for (UInt i = 0; i < spatial_dimension; ++i) { \ var[i] = var##_bounds[i]; \ } \ } \ DECLARE_GET_BOUND(Lower, lower) DECLARE_GET_BOUND(Upper, upper) DECLARE_GET_BOUND(LocalLower, local_lower) DECLARE_GET_BOUND(LocalUpper, local_upper) #undef DECLARE_GET_BOUND /* -------------------------------------------------------------------------- */ inline void Mesh::addConnectivityType(const ElementType & type){ getConnectivityPointer(type); } /* -------------------------------------------------------------------------- */ inline bool Mesh::isPureGhostNode(UInt n) const { return nodes_type ? (*nodes_type)(n) == -3 : false; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isLocalOrMasterNode(UInt n) const { return nodes_type ? (*nodes_type)(n) == -2 || (*nodes_type)(n) == -1 : true; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isLocalNode(UInt n) const { return nodes_type ? (*nodes_type)(n) == -1 : true; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isMasterNode(UInt n) const { return nodes_type ? (*nodes_type)(n) == -2 : false; } /* -------------------------------------------------------------------------- */ inline bool Mesh::isSlaveNode(UInt n) const { return nodes_type ? (*nodes_type)(n) >= 0 : false; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNodeGlobalId(UInt local_id) const { return nodes_global_ids ? (*nodes_global_ids)(local_id) : local_id; } /* -------------------------------------------------------------------------- */ inline UInt Mesh::getNbGlobalNodes() const { return nodes_global_ids ? nb_global_nodes : nodes->getSize(); } /* -------------------------------------------------------------------------- */ inline Int Mesh::getNodeType(UInt local_id) const { return nodes_type ? (*nodes_type)(local_id) : -1; } diff --git a/src/mesh_utils/mesh_io/mesh_io_msh.cc b/src/mesh_utils/mesh_io/mesh_io_msh.cc index 75083365d..143dbb0c3 100644 --- a/src/mesh_utils/mesh_io/mesh_io_msh.cc +++ b/src/mesh_utils/mesh_io/mesh_io_msh.cc @@ -1,473 +1,477 @@ /** * @file mesh_io_msh.cc * * @author Nicolas Richart * * @date Fri Jun 18 11:47:19 2010 * * @brief Read/Write for MSH files generated by gmsh * * @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 . * */ /* ----------------------------------------------------------------------------- Version (Legacy) 1.0 $NOD number-of-nodes node-number x-coord y-coord z-coord ... $ENDNOD $ELM number-of-elements elm-number elm-type reg-phys reg-elem number-of-nodes node-number-list ... $ENDELM ----------------------------------------------------------------------------- Version 2.1 $MeshFormat version-number file-type data-size $EndMeshFormat $Nodes number-of-nodes node-number x-coord y-coord z-coord ... $EndNodes $Elements number-of-elements elm-number elm-type number-of-tags < tag > ... node-number-list ... $EndElements $PhysicalNames number-of-names physical-dimension physical-number "physical-name" ... $EndPhysicalNames $NodeData number-of-string-tags < "string-tag" > ... number-of-real-tags < real-tag > ... number-of-integer-tags < integer-tag > ... node-number value ... ... $EndNodeData $ElementData number-of-string-tags < "string-tag" > ... number-of-real-tags < real-tag > ... number-of-integer-tags < integer-tag > ... elm-number value ... ... $EndElementData $ElementNodeData number-of-string-tags < "string-tag" > ... number-of-real-tags < real-tag > ... number-of-integer-tags < integer-tag > ... elm-number number-of-nodes-per-element value ... ... $ElementEndNodeData ----------------------------------------------------------------------------- elem-type 1: 2-node line. 2: 3-node triangle. 3: 4-node quadrangle. 4: 4-node tetrahedron. 5: 8-node hexahedron. 6: 6-node prism. 7: 5-node pyramid. 8: 3-node second order line 9: 6-node second order triangle 10: 9-node second order quadrangle 11: 10-node second order tetrahedron 12: 27-node second order hexahedron 13: 18-node second order prism 14: 14-node second order pyramid 15: 1-node point. 16: 8-node second order quadrangle 17: 20-node second order hexahedron 18: 15-node second order prism 19: 13-node second order pyramid 20: 9-node third order incomplete triangle 21: 10-node third order triangle 22: 12-node fourth order incomplete triangle 23: 15-node fourth order triangle 24: 15-node fifth order incomplete triangle 25: 21-node fifth order complete triangle 26: 4-node third order edge 27: 5-node fourth order edge 28: 6-node fifth order edge 29: 20-node third order tetrahedron 30: 35-node fourth order tetrahedron 31: 56-node fifth order tetrahedron -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include "mesh_io_msh.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* Methods Implentations */ /* -------------------------------------------------------------------------- */ MeshIOMSH::MeshIOMSH() { canReadSurface = true; canReadExtendedData = true; _msh_nodes_per_elem[_msh_not_defined ] = 0; _msh_nodes_per_elem[_msh_segment_2 ] = 2; _msh_nodes_per_elem[_msh_triangle_3 ] = 3; _msh_nodes_per_elem[_msh_quadrangle_4 ] = 4; _msh_nodes_per_elem[_msh_tetrahedron_4 ] = 4; _msh_nodes_per_elem[_msh_hexahedron_8 ] = 8; _msh_nodes_per_elem[_msh_prism_1 ] = 1; _msh_nodes_per_elem[_msh_pyramid_1 ] = 1; _msh_nodes_per_elem[_msh_segment_3 ] = 3; _msh_nodes_per_elem[_msh_triangle_6 ] = 6; _msh_nodes_per_elem[_msh_quadrangle_9 ] = 9; _msh_nodes_per_elem[_msh_tetrahedron_10] = 10; _msh_nodes_per_elem[_msh_hexahedron_27 ] = 27; _msh_nodes_per_elem[_msh_prism_18 ] = 18; _msh_nodes_per_elem[_msh_pyramid_14 ] = 14; _msh_nodes_per_elem[_msh_point ] = 1; _msh_nodes_per_elem[_msh_quadrangle_8 ] = 8; _msh_to_akantu_element_types[_msh_not_defined ] = _not_defined; _msh_to_akantu_element_types[_msh_segment_2 ] = _segment_2; _msh_to_akantu_element_types[_msh_triangle_3 ] = _triangle_3; _msh_to_akantu_element_types[_msh_quadrangle_4 ] = _quadrangle_4; _msh_to_akantu_element_types[_msh_tetrahedron_4 ] = _tetrahedron_4; _msh_to_akantu_element_types[_msh_hexahedron_8 ] = _hexahedron_8; _msh_to_akantu_element_types[_msh_prism_1 ] = _not_defined; _msh_to_akantu_element_types[_msh_pyramid_1 ] = _not_defined; _msh_to_akantu_element_types[_msh_segment_3 ] = _segment_3; _msh_to_akantu_element_types[_msh_triangle_6 ] = _triangle_6; _msh_to_akantu_element_types[_msh_quadrangle_9 ] = _not_defined; _msh_to_akantu_element_types[_msh_tetrahedron_10] = _tetrahedron_10; _msh_to_akantu_element_types[_msh_hexahedron_27 ] = _not_defined; _msh_to_akantu_element_types[_msh_prism_18 ] = _not_defined; _msh_to_akantu_element_types[_msh_pyramid_14 ] = _not_defined; _msh_to_akantu_element_types[_msh_point ] = _not_defined; _msh_to_akantu_element_types[_msh_quadrangle_8 ] = _quadrangle_8; _akantu_to_msh_element_types[_not_defined ] = _msh_not_defined; _akantu_to_msh_element_types[_segment_2 ] = _msh_segment_2; _akantu_to_msh_element_types[_segment_3 ] = _msh_segment_3; _akantu_to_msh_element_types[_triangle_3 ] = _msh_triangle_3; _akantu_to_msh_element_types[_triangle_6 ] = _msh_triangle_6; _akantu_to_msh_element_types[_tetrahedron_4 ] = _msh_tetrahedron_4; _akantu_to_msh_element_types[_tetrahedron_10 ] = _msh_tetrahedron_10; _akantu_to_msh_element_types[_quadrangle_4 ] = _msh_quadrangle_4; _akantu_to_msh_element_types[_quadrangle_8 ] = _msh_quadrangle_8; _akantu_to_msh_element_types[_hexahedron_8 ] = _msh_hexahedron_8; _akantu_to_msh_element_types[_point ] = _msh_point; _akantu_to_msh_element_types[_bernoulli_beam_2] = _msh_segment_2; std::map::iterator it; for(it = _akantu_to_msh_element_types.begin(); it != _akantu_to_msh_element_types.end(); ++it) { UInt nb_nodes = _msh_nodes_per_elem[it->second]; UInt * tmp = new UInt[nb_nodes]; for (UInt i = 0; i < nb_nodes; ++i) { tmp[i] = i; } switch(it->first) { case _tetrahedron_10: tmp[8] = 9; tmp[9] = 8; break; default: //nothing to change break; } _read_order[it->first] = tmp; } } /* -------------------------------------------------------------------------- */ MeshIOMSH::~MeshIOMSH() { std::map::iterator it; for(it = _akantu_to_msh_element_types.begin(); it != _akantu_to_msh_element_types.end(); ++it) { delete [] _read_order[it->first]; } } /* -------------------------------------------------------------------------- */ void MeshIOMSH::read(const std::string & filename, Mesh & mesh) { std::ifstream infile; infile.open(filename.c_str()); std::string line; UInt first_node_number = std::numeric_limits::max(), last_node_number = 0, file_format = 1, current_line = 0; if(!infile.good()) { AKANTU_DEBUG_ERROR("Cannot open file " << filename); } while(infile.good()) { std::getline(infile, line); current_line++; /// read the header if(line == "$MeshFormat") { std::getline(infile, line); /// the format line std::stringstream sstr(line); std::string version; sstr >> version; Int format; sstr >> format; if(format != 0) AKANTU_DEBUG_ERROR("This reader can only read ASCII files."); std::getline(infile, line); /// the end of block line current_line += 2; file_format = 2; } /// read all nodes if(line == "$Nodes" || line == "$NOD") { UInt nb_nodes; std::getline(infile, line); std::stringstream sstr(line); sstr >> nb_nodes; current_line++; Vector & nodes = const_cast &>(mesh.getNodes()); nodes.resize(nb_nodes); mesh.nb_global_nodes = nb_nodes; UInt index; Real coord[3]; UInt spatial_dimension = nodes.getNbComponent(); /// for each node, read the coordinates for(UInt i = 0; i < nb_nodes; ++i) { UInt offset = i * spatial_dimension; std::getline(infile, line); std::stringstream sstr_node(line); sstr_node >> index >> coord[0] >> coord[1] >> coord[2]; current_line++; first_node_number = std::min(first_node_number,index); last_node_number = std::max(last_node_number, index); /// read the coordinates for(UInt j = 0; j < spatial_dimension; ++j) nodes.values[offset + j] = coord[j]; } std::getline(infile, line); /// the end of block line } /// read all elements if(line == "$Elements" || line == "$ELM") { UInt nb_elements; UInt * read_order = NULL; std::getline(infile, line); std::stringstream sstr(line); sstr >> nb_elements; current_line++; Int index; UInt msh_type; ElementType akantu_type, akantu_type_old = _not_defined; Vector *connectivity = NULL; UInt node_per_element = 0; for(UInt i = 0; i < nb_elements; ++i) { std::getline(infile, line); std::stringstream sstr_elem(line); current_line++; sstr_elem >> index; sstr_elem >> msh_type; /// get the connectivity vector depending on the element type akantu_type = _msh_to_akantu_element_types[(MSHElementType) msh_type]; if(akantu_type == _not_defined) { AKANTU_DEBUG_WARNING("Unsuported element kind " << msh_type << " at line " << current_line); continue; } if(akantu_type != akantu_type_old) { connectivity = mesh.getConnectivityPointer(akantu_type); connectivity->resize(0); node_per_element = connectivity->getNbComponent(); akantu_type_old = akantu_type; read_order = _read_order[akantu_type]; } /// read tags informations if(file_format == 2) { UInt nb_tags; sstr_elem >> nb_tags; for(UInt j = 0; j < nb_tags; ++j) { Int tag; sstr_elem >> tag; std::stringstream sstr_tag_name; sstr_tag_name << "tag_" << j; Vector * data = mesh.getUIntDataPointer(akantu_type, sstr_tag_name.str(), _not_ghost); data->push_back(tag); } } else if (file_format == 1) { Int tag; sstr_elem >> tag; //reg-phys std::string tag_name = "tag_0"; Vector * data = mesh.getUIntDataPointer(akantu_type, tag_name, _not_ghost); data->push_back(tag); sstr_elem >> tag; //reg-elem tag_name = "tag_1"; data = mesh.getUIntDataPointer(akantu_type, tag_name, _not_ghost); data->push_back(tag); sstr_elem >> tag; //number-of-nodes } UInt local_connect[node_per_element]; for(UInt j = 0; j < node_per_element; ++j) { UInt node_index; sstr_elem >> node_index; AKANTU_DEBUG_ASSERT(node_index <= last_node_number, "Node number not in range : line " << current_line); node_index -= first_node_number; local_connect[read_order[j]] = node_index; } connectivity->push_back(local_connect); } std::getline(infile, line); /// the end of block line } if((line[0] == '$') && (line.find("End") == std::string::npos)) { AKANTU_DEBUG_WARNING("Unsuported block_kind " << line << " at line " << current_line); } } mesh.updateTypesOffsets(_not_ghost); infile.close(); } /* -------------------------------------------------------------------------- */ void MeshIOMSH::write(const std::string & filename, const Mesh & mesh) { std::ofstream outfile; const Vector & nodes = mesh.getNodes(); outfile.open(filename.c_str()); outfile << "$MeshFormat" << std::endl; outfile << "2.1 0 8" << std::endl;; outfile << "$EndMeshFormat" << std::endl;; outfile << std::setprecision(std::numeric_limits::digits10); outfile << "$Nodes" << std::endl;; outfile << nodes.getSize() << std::endl;; outfile << std::uppercase; for(UInt i = 0; i < nodes.getSize(); ++i) { Int offset = i * nodes.getNbComponent(); outfile << i+1; for(UInt j = 0; j < nodes.getNbComponent(); ++j) { outfile << " " << nodes.values[offset + j]; } if(nodes.getNbComponent() == 2) outfile << " " << 0.; outfile << std::endl;; } outfile << std::nouppercase; outfile << "$EndNodes" << std::endl;; outfile << "$Elements" << std::endl;; Mesh::type_iterator it = mesh.firstType(); Mesh::type_iterator end = mesh.lastType(); Int nb_elements = 0; for(; it != end; ++it) { const Vector & connectivity = mesh.getConnectivity(*it, _not_ghost); nb_elements += connectivity.getSize(); } outfile << nb_elements << std::endl; UInt element_idx = 1; for(it = mesh.firstType(); it != end; ++it) { ElementType type = *it; const Vector & connectivity = mesh.getConnectivity(type, _not_ghost); + const Vector & data_tag_0 = mesh.getUIntData(type, "tag_0", _not_ghost); + const Vector & data_tag_1 = mesh.getUIntData(type, "tag_1", _not_ghost); for(UInt i = 0; i < connectivity.getSize(); ++i) { UInt offset = i * connectivity.getNbComponent(); - outfile << element_idx << " " << _akantu_to_msh_element_types[type] << " 3 1 1 0"; /// \todo write the real data in the file + outfile << element_idx << " " << _akantu_to_msh_element_types[type] + << " 3 " << data_tag_0(i) << " " << data_tag_1(i) + << " 0"; /// \todo write the real data in the file for(UInt j = 0; j < connectivity.getNbComponent(); ++j) { outfile << " " << connectivity.values[offset + j] + 1; } outfile << std::endl; element_idx++; } } outfile << "$EndElements" << std::endl;; outfile.close(); } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/mesh_utils/mesh_partition.cc b/src/mesh_utils/mesh_partition.cc index cdd3f83f8..0cf8081bb 100644 --- a/src/mesh_utils/mesh_partition.cc +++ b/src/mesh_utils/mesh_partition.cc @@ -1,304 +1,295 @@ /** * @file mesh_partition.cc * * @author David Simon Kammer * @author Nicolas Richart * * @date Tue Aug 17 16:19:43 2010 * * @brief implementation of common part of all partitioner * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "mesh_partition.hh" #include "mesh_utils.hh" #include "aka_types.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ MeshPartition::MeshPartition(const Mesh & mesh, UInt spatial_dimension, const MemoryID & memory_id) : Memory(memory_id), id("MeshPartitioner"), mesh(mesh), spatial_dimension(spatial_dimension), - partitions ("partition" , id), - ghost_partitions ("ghost_partition" , id), - ghost_partitions_offset("ghost_partition_offset", id) { + partitions ("partition" , id, memory_id), + ghost_partitions ("ghost_partition" , id, memory_id), + ghost_partitions_offset("ghost_partition_offset", id, memory_id) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MeshPartition::~MeshPartition() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** - * conversion in c++ of a the GENDUALMETIS (mesh.c) function wrote by George in + * conversion in c++ of the GENDUALMETIS (mesh.c) function wrote by George in * Metis (University of Minnesota) */ -void MeshPartition::buildDualGraph(Vector & dxadj, Vector & dadjncy) { +void MeshPartition::buildDualGraph(Vector & dxadj, Vector & dadjncy, + Vector & edge_loads, + const EdgeLoadFunctor & edge_load_func, + const Vector & pairs) { AKANTU_DEBUG_IN(); + // tweak mesh; const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); UInt nb_types = type_list.size(); UInt nb_good_types = 0; UInt nb_nodes_per_element[nb_types]; UInt nb_nodes_per_element_p1[nb_types]; UInt magic_number[nb_types]; - UInt * conn_val[nb_types]; + // UInt * conn_val[nb_types]; UInt nb_element[nb_types]; + Vector * conn[nb_types]; + Vector * conn_tmp[nb_types]; + + Vector lin_to_element; + + Element el; + el.ghost_type = _not_ghost; + Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; if(Mesh::getSpatialDimension(type) != mesh.getSpatialDimension()) continue; + el.type = type; ElementType type_p1 = Mesh::getP1ElementType(type); nb_nodes_per_element[nb_good_types] = Mesh::getNbNodesPerElement(type); nb_nodes_per_element_p1[nb_good_types] = Mesh::getNbNodesPerElement(type_p1); - conn_val[nb_good_types] = mesh.getConnectivity(type, _not_ghost).values; nb_element[nb_good_types] = mesh.getConnectivity(type, _not_ghost).getSize(); magic_number[nb_good_types] = Mesh::getNbNodesPerElement(Mesh::getFacetElementType(type_p1)); - nb_good_types++; - } + conn[nb_good_types] = &const_cast &>(mesh.getConnectivity(type, _not_ghost)); - // Vector node_offset; - // Vector node_index; + for (UInt i = 0; i < nb_element[nb_good_types]; ++i) { + el.element = i; + lin_to_element.push_back(el); + } - CSR node_to_elem; - // MeshUtils::buildNode2Elements(mesh, node_offset, node_index); - MeshUtils::buildNode2Elements(mesh, node_to_elem); + if(pairs.getSize() != 0) { + conn_tmp[nb_good_types] = new Vector(mesh.getConnectivity(type, _not_ghost)); + for (UInt i = 0; i < pairs.getSize(); ++i) { + for (UInt el = 0; el < nb_element[nb_good_types]; ++el) { + for (UInt n = 0; n < nb_nodes_per_element[nb_good_types]; ++n) { + if(pairs(i, 1) == (*conn[nb_good_types])(el, n)) + (*conn[nb_good_types])(el, n) = pairs(i, 0); + } + } + } + } - // UInt * node_offset_val = node_offset.values; - // UInt * node_index_val = node_index.values; + nb_good_types++; + } + CSR node_to_elem; + + MeshUtils::buildNode2Elements(mesh, node_to_elem); UInt nb_total_element = 0; UInt nb_total_node_element = 0; for (UInt t = 0; t < nb_good_types; ++t) { nb_total_element += nb_element[t]; nb_total_node_element += nb_element[t]*nb_nodes_per_element_p1[t]; } dxadj.resize(nb_total_element + 1); - Int * dxadj_val = dxadj.values; /// initialize the dxadj array for (UInt t = 0, linerized_el = 0; t < nb_good_types; ++t) for (UInt el = 0; el < nb_element[t]; ++el, ++linerized_el) - dxadj_val[linerized_el] = nb_nodes_per_element_p1[t]; - + dxadj(linerized_el) = nb_nodes_per_element_p1[t]; /// convert the dxadj_val array in a csr one - for (UInt i = 1; i < nb_total_element; ++i) dxadj_val[i] += dxadj_val[i-1]; - for (UInt i = nb_total_element; i > 0; --i) dxadj_val[i] = dxadj_val[i-1]; - dxadj_val[0] = 0; + for (UInt i = 1; i < nb_total_element; ++i) dxadj(i) += dxadj(i-1); + for (UInt i = nb_total_element; i > 0; --i) dxadj(i) = dxadj(i-1); + dxadj(0) = 0; - dadjncy.resize(2*dxadj_val[nb_total_element]); - Int * dadjncy_val = dadjncy.values; + dadjncy.resize(2*dxadj(nb_total_element)); + edge_loads.resize(2*dxadj(nb_total_element)); /// weight map to determine adjacency unordered_map::type weight_map; - // UInt index[200], weight[200]; /// key, value - // UInt mask = (1 << 11) - 1; /// hash function - // Int * mark = new Int[mask + 1]; /// collision detector - // for (UInt i = 0; i < mask + 1; ++i) mark[i] = -1; - - for (UInt t = 0, linerized_el = 0; t < nb_good_types; ++t) { for (UInt el = 0; el < nb_element[t]; ++el, ++linerized_el) { - UInt el_offset = el*nb_nodes_per_element[t]; - /// fill the weight map - // UInt m = 0; for (UInt n = 0; n < nb_nodes_per_element_p1[t]; ++n) { - UInt node = conn_val[t][el_offset + n]; + UInt node = (*conn[t])(el, n); CSR::iterator k; for (k = node_to_elem.rbegin(node); k != node_to_elem.rend(node); --k) { - // for (UInt k = node_offset_val[node + 1] - 1; - // k >= node_offset_val[node]; - // --k) { - UInt current_el = *k;//node_index_val[k]; + UInt current_el = *k; if(current_el <= linerized_el) break; unordered_map::type::iterator it_w; it_w = weight_map.find(current_el); if(it_w == weight_map.end()) { weight_map[current_el] = 1; } else { it_w->second++; } - - // UInt mark_offset = current_el & mask; - // Int current_mark = mark[mark_offset]; - // if(current_mark == -1) { /// if element not in map - // index[m] = current_el; - // weight[m] = 1; - // mark[mark_offset] = m++; - // } else if (index[current_mark] == current_el) { /// if element already in map - // weight[current_mark]++; - // } else { /// if element already in map but collision in the keys - // UInt i; - // for (i = 0; i < m; ++i) { - // if(index[i] == current_el) { - // weight[i]++; - // break; - // } - // } - // if(i == m) { - // index[m] = current_el; - // weight[m++] = 1; - // } - // } } } /// each element with a weight of the size of a facet are adjacent unordered_map::type::iterator it_w; for(it_w = weight_map.begin(); it_w != weight_map.end(); ++it_w) { if(it_w->second == magic_number[t]) { UInt adjacent_el = it_w->first; - UInt index_adj = dxadj_val[adjacent_el ]++; - UInt index_lin = dxadj_val[linerized_el]++; + UInt index_adj = dxadj(adjacent_el )++; + UInt index_lin = dxadj(linerized_el)++; - dadjncy_val[index_lin] = adjacent_el; - dadjncy_val[index_adj] = linerized_el; + dadjncy(index_lin) = adjacent_el; + dadjncy(index_adj) = linerized_el; } } - // for (UInt n = 0; n < m; ++n) { - // if(weight[n] == magic_number[t]) { - // UInt adjacent_el = index[n]; - // dadjncy_val[dxadj_val[linerized_el]++] = adjacent_el; - // dadjncy_val[dxadj_val[adjacent_el ]++] = linerized_el; - // } - // mark[index[n] & mask] = -1; - // } - weight_map.clear(); } } + if(pairs.getSize() != 0) { + for (UInt i = 0; i < nb_good_types; ++i) { + conn[i]->copy(*conn_tmp[i]); + delete conn_tmp[i]; + } + } + Int k_start = 0; for (UInt t = 0, linerized_el = 0, j = 0; t < nb_good_types; ++t) for (UInt el = 0; el < nb_element[t]; ++el, ++linerized_el) { - for (Int k = k_start; k < dxadj_val[linerized_el]; ++k, ++j) - dadjncy_val[j] = dadjncy_val[k]; - dxadj_val[linerized_el] = j; + for (Int k = k_start; k < dxadj(linerized_el); ++k, ++j) + dadjncy(j) = dadjncy(k); + dxadj(linerized_el) = j; k_start += nb_nodes_per_element_p1[t]; } - for (UInt i = nb_total_element; i > 0; --i) dxadj_val[i] = dxadj_val[i-1]; - dxadj_val[0] = 0; + for (UInt i = nb_total_element; i > 0; --i) dxadj(i) = dxadj(i - 1); + dxadj(0) = 0; - // delete [] mark; + UInt adj = 0; + for (UInt i = 0; i < nb_total_element; ++i) { + UInt nb_adj = dxadj(i + 1) - dxadj(i); + for (UInt j = 0; j < nb_adj; ++j, ++adj) { + Int el_adj_id = dadjncy(dxadj(i) + j); + Element el = lin_to_element(i); + Element el_adj = lin_to_element(el_adj_id); + + Int load = edge_load_func(el, el_adj); + edge_loads(adj) = load; + } + } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshPartition::fillPartitionInformation(const Mesh & mesh, const Int * linearized_partitions) { AKANTU_DEBUG_IN(); - // Vector node_offset; - // Vector node_index; CSR node_to_elem; - // MeshUtils::buildNode2Elements(mesh, node_offset, node_index); MeshUtils::buildNode2Elements(mesh, node_to_elem); - // UInt * node_offset_val = node_offset.values; - // UInt * node_index_val = node_index.values; Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator end = mesh.lastType(spatial_dimension); UInt linearized_el = 0; for(; it != end; ++it) { ElementType type = *it; UInt nb_element = mesh.getNbElement(*it); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); partitions .alloc(nb_element, 1, type, _not_ghost); ghost_partitions_offset.alloc(nb_element + 1, 1, type, _ghost); ghost_partitions .alloc(0, 1, type, _ghost); const Vector & connectivity = mesh.getConnectivity(type, _not_ghost); for (UInt el = 0; el < nb_element; ++el, ++linearized_el) { UInt part = linearized_partitions[linearized_el]; partitions(type, _not_ghost)(el) = part; std::list list_adj_part; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = connectivity.values[el * nb_nodes_per_element + n]; CSR::iterator ne; - // for (UInt ne = node_offset_val[node]; ne < node_offset_val[node + 1]; ++ne) { for (ne = node_to_elem.begin(node); ne != node_to_elem.end(node); ++ne) { - UInt adj_el = *ne;//node_index_val[ne]; + UInt adj_el = *ne; UInt adj_part = linearized_partitions[adj_el]; if(part != adj_part) { list_adj_part.push_back(adj_part); } } } list_adj_part.sort(); list_adj_part.unique(); for(std::list::iterator adj_it = list_adj_part.begin(); adj_it != list_adj_part.end(); ++adj_it) { ghost_partitions(type, _ghost).push_back(*adj_it); ghost_partitions_offset(type, _ghost)(el)++; } } /// convert the ghost_partitions_offset array in an offset array Vector & ghost_partitions_offset_ptr = ghost_partitions_offset(type, _ghost); for (UInt i = 1; i < nb_element; ++i) ghost_partitions_offset_ptr(i) += ghost_partitions_offset_ptr(i-1); for (UInt i = nb_element; i > 0; --i) ghost_partitions_offset_ptr(i) = ghost_partitions_offset_ptr(i-1); ghost_partitions_offset_ptr(0) = 0; } AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/src/mesh_utils/mesh_partition.hh b/src/mesh_utils/mesh_partition.hh index 592700c67..6cd862c36 100644 --- a/src/mesh_utils/mesh_partition.hh +++ b/src/mesh_utils/mesh_partition.hh @@ -1,132 +1,153 @@ /** * @file mesh_partition.hh * * @author David Simon Kammer * @author Nicolas Richart * * @date Mon Aug 16 13:17:22 2010 * * @brief tools to partitionate a mesh * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_PARTITION_HH__ #define __AKANTU_MESH_PARTITION_HH__ /* -------------------------------------------------------------------------- */ #include "aka_memory.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ class MeshPartition : protected Memory { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MeshPartition(const Mesh & mesh, UInt spatial_dimension, const MemoryID & memory_id = 0); virtual ~MeshPartition(); + class EdgeLoadFunctor { + public: + virtual inline Int operator()(__attribute__((unused)) const Element & el1, + __attribute__((unused)) const Element & el2) const = 0; + }; + + class ConstEdgeLoadFunctor : public EdgeLoadFunctor { + public: + virtual inline Int operator()(__attribute__((unused)) const Element & el1, + __attribute__((unused)) const Element & el2) const { + return 1; + } + }; + + /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: - virtual void partitionate(UInt nb_part) = 0; + virtual void partitionate(UInt nb_part, + const EdgeLoadFunctor & edge_load_func = ConstEdgeLoadFunctor(), + const Vector & pairs = Vector()) = 0; virtual void reorder() = 0; /// function to print the contain of the class //virtual void printself(std::ostream & stream, int indent = 0) const = 0; /// fill the partitions array with a given linearized partition information void fillPartitionInformation(const Mesh & mesh, const Int * linearized_partitions); protected: /// build the dual graph of the mesh, for all element of spatial_dimension - void buildDualGraph(Vector & dxadj, Vector & dadjncy); + void buildDualGraph(Vector & dxadj, Vector & dadjncy, + Vector & edge_loads, + const EdgeLoadFunctor & edge_load_func, + const Vector & pairs); + /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Partitions, partitions, const ByElementTypeUInt &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Partition, partitions, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(GhostPartition, ghost_partitions, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(GhostPartitionOffset, ghost_partitions_offset, UInt); AKANTU_GET_MACRO(NbPartition, nb_partitions, UInt); AKANTU_SET_MACRO(NbPartition, nb_partitions, UInt); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id std::string id; /// the mesh to partition const Mesh & mesh; /// dimension of the elements to consider in the mesh UInt spatial_dimension; /// number of partitions UInt nb_partitions; /// partition numbers ByElementTypeUInt partitions; ByElementTypeUInt ghost_partitions; ByElementTypeUInt ghost_partitions_offset; Vector * permutation; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "mesh_partition_inline_impl.cc" /// standard output stream operator // inline std::ostream & operator <<(std::ostream & stream, const MeshPartition & _this) // { // _this.printself(stream); // return stream; // } __END_AKANTU__ #ifdef AKANTU_USE_SCOTCH #include "mesh_partition_scotch.hh" #endif #endif /* __AKANTU_MESH_PARTITION_HH__ */ diff --git a/src/mesh_utils/mesh_partition/mesh_partition_scotch.cc b/src/mesh_utils/mesh_partition/mesh_partition_scotch.cc index e95b20e49..47a52b818 100644 --- a/src/mesh_utils/mesh_partition/mesh_partition_scotch.cc +++ b/src/mesh_utils/mesh_partition/mesh_partition_scotch.cc @@ -1,459 +1,461 @@ /** * @file mesh_partition_scotch.cc * * @author David Simon Kammer * @author Nicolas Richart * * @date Tue Aug 17 16:19:43 2010 * * @brief implementation of the MeshPartitionScotch class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include "mesh_partition_scotch.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ MeshPartitionScotch::MeshPartitionScotch(const Mesh & mesh, UInt spatial_dimension, const MemoryID & memory_id) : MeshPartition(mesh, spatial_dimension, memory_id) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SCOTCH_Mesh * MeshPartitionScotch::createMesh() { AKANTU_DEBUG_IN(); UInt spatial_dimension = mesh.getSpatialDimension(); UInt nb_nodes = mesh.getNbNodes(); UInt total_nb_element = 0; UInt nb_edge = 0; Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator end = mesh.lastType(spatial_dimension); for(; it != end; ++it) { ElementType type = *it; UInt nb_element = mesh.getNbElement(type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); total_nb_element += nb_element; nb_edge += nb_element * nb_nodes_per_element; } SCOTCH_Num vnodbas = 0; SCOTCH_Num vnodnbr = nb_nodes; SCOTCH_Num velmbas = vnodnbr; SCOTCH_Num velmnbr = total_nb_element; SCOTCH_Num * verttab = new SCOTCH_Num[vnodnbr + velmnbr + 1]; SCOTCH_Num * vendtab = verttab + 1; SCOTCH_Num * velotab = NULL; SCOTCH_Num * vnlotab = NULL; SCOTCH_Num * vlbltab = NULL; memset(verttab, 0, (vnodnbr + velmnbr + 1) * sizeof(SCOTCH_Num)); it = mesh.firstType(spatial_dimension); for(; it != end; ++it) { ElementType type = *it; if(Mesh::getSpatialDimension(type) != spatial_dimension) continue; UInt nb_element = mesh.getNbElement(type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const Vector & connectivity = mesh.getConnectivity(type, _not_ghost); /// count number of occurrence of each node for (UInt el = 0; el < nb_element; ++el) { UInt * conn_val = connectivity.values + el * nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) { verttab[*(conn_val++)]++; } } } /// convert the occurrence array in a csr one for (UInt i = 1; i < nb_nodes; ++i) verttab[i] += verttab[i-1]; for (UInt i = nb_nodes; i > 0; --i) verttab[i] = verttab[i-1]; verttab[0] = 0; /// rearrange element to get the node-element list SCOTCH_Num edgenbr = verttab[vnodnbr] + nb_edge; SCOTCH_Num * edgetab = new SCOTCH_Num[edgenbr]; UInt linearized_el = 0; it = mesh.firstType(spatial_dimension); for(; it != end; ++it) { ElementType type = *it; UInt nb_element = mesh.getNbElement(type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const Vector & connectivity = mesh.getConnectivity(type, _not_ghost); for (UInt el = 0; el < nb_element; ++el, ++linearized_el) { UInt * conn_val = connectivity.values + el * nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) edgetab[verttab[*(conn_val++)]++] = linearized_el + velmbas; } } for (UInt i = nb_nodes; i > 0; --i) verttab[i] = verttab[i-1]; verttab[0] = 0; SCOTCH_Num * verttab_tmp = verttab + vnodnbr + 1; SCOTCH_Num * edgetab_tmp = edgetab + verttab[vnodnbr]; it = mesh.firstType(spatial_dimension); for(; it != end; ++it) { ElementType type = *it; UInt nb_element = mesh.getNbElement(type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const Vector & connectivity = mesh.getConnectivity(type, _not_ghost); for (UInt el = 0; el < nb_element; ++el) { *verttab_tmp = *(verttab_tmp - 1) + nb_nodes_per_element; verttab_tmp++; UInt * conn = connectivity.values + el * nb_nodes_per_element; for (UInt i = 0; i < nb_nodes_per_element; ++i) { *(edgetab_tmp++) = *(conn++) + vnodbas; } } } SCOTCH_Mesh * meshptr = new SCOTCH_Mesh; SCOTCH_meshInit(meshptr); SCOTCH_meshBuild(meshptr, velmbas, vnodbas, velmnbr, vnodnbr, verttab, vendtab, velotab, vnlotab, vlbltab, edgenbr, edgetab); /// Check the mesh AKANTU_DEBUG_ASSERT(SCOTCH_meshCheck(meshptr) == 0, "Scotch mesh is not consistent"); #ifndef AKANTU_NDEBUG if (AKANTU_DEBUG_TEST(dblDump)) { /// save initial graph FILE *fmesh = fopen("ScotchMesh.msh", "w"); SCOTCH_meshSave(meshptr, fmesh); fclose(fmesh); /// write geometry file std::ofstream fgeominit; fgeominit.open("ScotchMesh.xyz"); fgeominit << spatial_dimension << std::endl << nb_nodes << std::endl; const Vector & nodes = mesh.getNodes(); Real * nodes_val = nodes.values; for (UInt i = 0; i < nb_nodes; ++i) { fgeominit << i << " "; for (UInt s = 0; s < spatial_dimension; ++s) fgeominit << *(nodes_val++) << " "; fgeominit << std::endl;; } fgeominit.close(); } #endif AKANTU_DEBUG_OUT(); return meshptr; } /* -------------------------------------------------------------------------- */ void MeshPartitionScotch::destroyMesh(SCOTCH_Mesh * meshptr) { AKANTU_DEBUG_IN(); SCOTCH_Num velmbas, vnodbas, vnodnbr, velmnbr, *verttab, *vendtab, *velotab, *vnlotab, *vlbltab, edgenbr, *edgetab, degrptr; SCOTCH_meshData(meshptr, &velmbas, &vnodbas, &velmnbr, &vnodnbr, &verttab, &vendtab, &velotab, &vnlotab, &vlbltab, &edgenbr, &edgetab, °rptr); delete [] verttab; delete [] edgetab; SCOTCH_meshExit(meshptr); delete meshptr; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -void MeshPartitionScotch::partitionate(UInt nb_part) { +void MeshPartitionScotch::partitionate(UInt nb_part, + const EdgeLoadFunctor & edge_load_func, + const Vector & pairs) { AKANTU_DEBUG_IN(); nb_partitions = nb_part; AKANTU_DEBUG_INFO("Partitioning the mesh " << mesh.getID() << " in " << nb_part << " parts."); Vector dxadj; Vector dadjncy; - - buildDualGraph(dxadj, dadjncy); + Vector edge_loads; + buildDualGraph(dxadj, dadjncy, edge_loads, edge_load_func, pairs); /// variables that will hold our structures in scotch format SCOTCH_Graph scotch_graph; SCOTCH_Strat scotch_strat; /// description number and arrays for struct mesh for scotch SCOTCH_Num baseval = 0; //base numbering for element and //nodes (0 -> C , 1 -> fortran) SCOTCH_Num vertnbr = dxadj.getSize() - 1; //number of vertexes SCOTCH_Num *parttab; //array of partitions - SCOTCH_Num edgenbr = dxadj.values[vertnbr]; //twice the number of "edges" + SCOTCH_Num edgenbr = dxadj(vertnbr); //twice the number of "edges" //(an "edge" bounds two nodes) - SCOTCH_Num * verttab = dxadj.values; //array of start indices in edgetab + SCOTCH_Num * verttab = dxadj.storage(); //array of start indices in edgetab SCOTCH_Num * vendtab = NULL; //array of after-last indices in edgetab SCOTCH_Num * velotab = NULL; //integer load associated with //every vertex ( optional ) - SCOTCH_Num * edlotab = NULL; //integer load associated with + SCOTCH_Num * edlotab = edge_loads.storage(); //integer load associated with //every edge ( optional ) - SCOTCH_Num * edgetab = dadjncy.values; // adjacency array of every vertex + SCOTCH_Num * edgetab = dadjncy.storage(); // adjacency array of every vertex SCOTCH_Num * vlbltab = NULL; // vertex label array (optional) /// Allocate space for Scotch arrays parttab = new SCOTCH_Num[vertnbr]; /// Initialize the strategy structure SCOTCH_stratInit (&scotch_strat); /// Initialize the graph structure SCOTCH_graphInit(&scotch_graph); /// Build the graph from the adjacency arrays SCOTCH_graphBuild(&scotch_graph, baseval, vertnbr, verttab, vendtab, velotab, vlbltab, edgenbr, edgetab, edlotab); #ifndef AKANTU_NDEBUG if (AKANTU_DEBUG_TEST(dblDump)) { /// save initial graph FILE *fgraphinit = fopen("GraphIniFile.grf", "w"); - SCOTCH_graphSave(&scotch_graph,fgraphinit); + SCOTCH_graphSave(&scotch_graph, fgraphinit); fclose(fgraphinit); /// write geometry file std::ofstream fgeominit; fgeominit.open("GeomIniFile.xyz"); fgeominit << spatial_dimension << std::endl << vertnbr << std::endl; const Vector & nodes = mesh.getNodes(); Mesh::type_iterator f_it = mesh.firstType(spatial_dimension); Mesh::type_iterator f_end = mesh.lastType(spatial_dimension); UInt out_linerized_el = 0; for(; f_it != f_end; ++f_it) { ElementType type = *f_it; UInt nb_element = mesh.getNbElement(*f_it); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const Vector & connectivity = mesh.getConnectivity(type); Real mid[spatial_dimension] ; for (UInt el = 0; el < nb_element; ++el) { memset(mid, 0, spatial_dimension*sizeof(Real)); for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = connectivity.values[nb_nodes_per_element * el + n]; for (UInt s = 0; s < spatial_dimension; ++s) mid[s] += ((Real) nodes.values[node * spatial_dimension + s]) / ((Real) nb_nodes_per_element); } fgeominit << out_linerized_el++ << " "; for (UInt s = 0; s < spatial_dimension; ++s) fgeominit << mid[s] << " "; fgeominit << std::endl;; } } fgeominit.close(); } #endif /// Check the graph AKANTU_DEBUG_ASSERT(SCOTCH_graphCheck(&scotch_graph) == 0, "Graph to partition is not consistent"); /// Partition the mesh SCOTCH_graphPart(&scotch_graph, nb_part, &scotch_strat, parttab); /// Check the graph AKANTU_DEBUG_ASSERT(SCOTCH_graphCheck(&scotch_graph) == 0, "Partitioned graph is not consistent"); #ifndef AKANTU_NDEBUG if (AKANTU_DEBUG_TEST(dblDump)) { /// save the partitioned graph FILE *fgraph=fopen("GraphFile.grf", "w"); SCOTCH_graphSave(&scotch_graph, fgraph); fclose(fgraph); /// save the partition map std::ofstream fmap; fmap.open("MapFile.map"); fmap << vertnbr << std::endl; for (Int i = 0; i < vertnbr; i++) fmap << i << " " << parttab[i] << std::endl; fmap.close(); } #endif /// free the scotch data structures SCOTCH_stratExit(&scotch_strat); SCOTCH_graphFree(&scotch_graph); SCOTCH_graphExit(&scotch_graph); fillPartitionInformation(mesh, parttab); delete [] parttab; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MeshPartitionScotch::reorder() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Reordering the mesh " << mesh.getID()); SCOTCH_Mesh * scotch_mesh = createMesh(); UInt nb_nodes = mesh.getNbNodes(); SCOTCH_Strat scotch_strat; // SCOTCH_Ordering scotch_order; SCOTCH_Num * permtab = new SCOTCH_Num[nb_nodes]; SCOTCH_Num * peritab = NULL; SCOTCH_Num cblknbr = 0; SCOTCH_Num * rangtab = NULL; SCOTCH_Num * treetab = NULL; /// Initialize the strategy structure SCOTCH_stratInit (&scotch_strat); SCOTCH_Graph scotch_graph; SCOTCH_graphInit(&scotch_graph); SCOTCH_meshGraph(scotch_mesh, &scotch_graph); #ifndef AKANTU_NDEBUG if (AKANTU_DEBUG_TEST(dblDump)) { FILE *fgraphinit = fopen("ScotchMesh.grf", "w"); SCOTCH_graphSave(&scotch_graph,fgraphinit); fclose(fgraphinit); } #endif /// Check the graph // AKANTU_DEBUG_ASSERT(SCOTCH_graphCheck(&scotch_graph) == 0, // "Mesh to Graph is not consistent"); SCOTCH_graphOrder(&scotch_graph, &scotch_strat, permtab, peritab, &cblknbr, rangtab, treetab); SCOTCH_graphExit(&scotch_graph); SCOTCH_stratExit(&scotch_strat); destroyMesh(scotch_mesh); /// Renumbering UInt spatial_dimension = mesh.getSpatialDimension(); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; Mesh::type_iterator it = mesh.firstType(0, gt); Mesh::type_iterator end = mesh.lastType(0, gt); for(; it != end; ++it) { ElementType type = *it; UInt nb_element = mesh.getNbElement(type, gt); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); const Vector & connectivity = mesh.getConnectivity(type, gt); UInt * conn = connectivity.values; for (UInt el = 0; el < nb_element * nb_nodes_per_element; ++el, ++conn) { *conn = permtab[*conn]; } } } /// \todo think of a in-place way to do it Real * new_coordinates = new Real[spatial_dimension * nb_nodes]; Real * old_coordinates = mesh.getNodes().values; for (UInt i = 0; i < nb_nodes; ++i) { memcpy(new_coordinates + permtab[i] * spatial_dimension, old_coordinates + i * spatial_dimension, spatial_dimension * sizeof(Real)); } memcpy(old_coordinates, new_coordinates, nb_nodes * spatial_dimension * sizeof(Real)); delete [] new_coordinates; delete [] permtab; AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/src/mesh_utils/mesh_partition/mesh_partition_scotch.hh b/src/mesh_utils/mesh_partition/mesh_partition_scotch.hh index 03716621f..c9dd331ef 100644 --- a/src/mesh_utils/mesh_partition/mesh_partition_scotch.hh +++ b/src/mesh_utils/mesh_partition/mesh_partition_scotch.hh @@ -1,110 +1,112 @@ /** * @file mesh_partition_scotch.hh * * @author David Simon Kammer * @author Nicolas Richart * * @date Mon Aug 16 13:17:22 2010 * * @brief mesh partitioning based on libScotch * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_PARTITION_SCOTCH_HH__ #define __AKANTU_MESH_PARTITION_SCOTCH_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh_partition.hh" #if ! defined(AKANTU_USE_PTSCOTCH) # ifndef AKANTU_SCOTCH_NO_EXTERN extern "C" { # endif //AKANTU_SCOTCH_NO_EXTERN # include # ifndef AKANTU_SCOTCH_NO_EXTERN } # endif //AKANTU_SCOTCH_NO_EXTERN #else //AKANTU_USE_PTSCOTCH # include #endif //AKANTU_USE_PTSCOTCH /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ class MeshPartitionScotch : public MeshPartition { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MeshPartitionScotch(const Mesh & mesh, UInt spatial_dimension, const MemoryID & memory_id = 0); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: - virtual void partitionate(UInt nb_part); + virtual void partitionate(UInt nb_part, + const EdgeLoadFunctor & edge_load_func = ConstEdgeLoadFunctor(), + const Vector & pairs = Vector()); virtual void reorder(); /// function to print the contain of the class //virtual void printself(std::ostream & stream, int indent = 0) const; private: SCOTCH_Mesh * createMesh(); void destroyMesh(SCOTCH_Mesh * meshptr); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ // #include "mesh_partition_scotch_inline_impl.cc" /// standard output stream operator // inline std::ostream & operator <<(std::ostream & stream, const MeshPartitionScotch & _this) // { // _this.printself(stream); // return stream; // } __END_AKANTU__ #endif /* __AKANTU_MESH_PARTITION_SCOTCH_HH__ */ diff --git a/src/model/solid_mechanics/contact/contact_rigid.cc b/src/model/solid_mechanics/contact/contact_rigid.cc index 1966a4925..f2d2db933 100644 --- a/src/model/solid_mechanics/contact/contact_rigid.cc +++ b/src/model/solid_mechanics/contact/contact_rigid.cc @@ -1,1223 +1,1241 @@ /** * @file contact_rigid.cc * * @author David Simon Kammer * * @date Wed Jan 19 14:45:41 2011 * * @brief Specialization of the contact structure for 3d contact in explicit * time scheme * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "contact_rigid.hh" #include "contact_search.hh" /* -------------------------------------------------------------------------- */ #include __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ ContactRigid::ImpactorInformationPerMaster::ImpactorInformationPerMaster(const Surface master_id, const UInt spatial_dimension) : master_id(master_id) { AKANTU_DEBUG_IN(); this->impactor_surfaces = new std::vector(0); this->active_impactor_nodes = new Vector(0,1); //this->master_element_offset = new Vector(0,1); this->master_element_type = new std::vector(0); this->master_normals = new Vector(0, spatial_dimension); this->node_is_sticking = new Vector(0,2); this->friction_forces = new Vector(0,spatial_dimension); this->stick_positions = new Vector(0,spatial_dimension); this->residual_forces = new Vector(0,spatial_dimension); this->previous_velocities = new Vector(0,spatial_dimension); this->friction_resistances = new Vector(0); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ ContactRigid::ImpactorInformationPerMaster::~ImpactorInformationPerMaster() { AKANTU_DEBUG_IN(); delete this->impactor_surfaces; delete this->active_impactor_nodes; // delete this->master_element_offset; delete this->master_element_type; delete this->master_normals; delete this->node_is_sticking; delete this->friction_forces; delete this->stick_positions; delete this->residual_forces; delete this->previous_velocities; delete this->friction_resistances; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ ContactRigid::ContactRigid(const SolidMechanicsModel & model, const ContactType & type, const ContactID & id, const MemoryID & memory_id) : - Contact(model, type, id, memory_id), spatial_dimension(model.getSpatialDimension()), mesh(model.getFEM().getMesh()), prakash(false), ref_velocity(0.), characterstic_length(0.) { + Contact(model, type, id, memory_id), spatial_dimension(model.getSpatialDimension()), mesh(model.getFEM().getMesh()), prakash(false), ref_velocity(0.), characterstic_length(0.), dedontney(false), t_star(0.) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ ContactRigid::~ContactRigid() { AKANTU_DEBUG_IN(); ContactRigid::SurfaceToImpactInfoMap::iterator it_imp; for (it_imp = this->impactors_information.begin(); it_imp != this->impactors_information.end(); ++it_imp) { delete it_imp->second; } this->impactors_information.clear(); this->friction_coefficient.clear(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactRigid::addMasterSurface(const Surface & master_surface) { AKANTU_DEBUG_IN(); Contact::addMasterSurface(master_surface); ImpactorInformationPerMaster * impactor_info = new ImpactorInformationPerMaster(master_surface, this->spatial_dimension); this->impactors_information[master_surface] = impactor_info; //this->friction_coefficient[master_surface] = NULL; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactRigid::addImpactorSurfaceToMasterSurface(const Surface & impactor_surface, const Surface & master_surface) { AKANTU_DEBUG_IN(); ContactRigid::SurfaceToImpactInfoMap::iterator it_imp; it_imp = this->impactors_information.find(master_surface); AKANTU_DEBUG_ASSERT(it_imp != this->impactors_information.end(), "The master surface: " << master_surface << "couldn't be found and erased in impactors_information map"); it_imp->second->impactor_surfaces->push_back(impactor_surface); /* for (UInt m=0; m < this->impactors_information.size(); ++m) { if (this->impactors_information.at(m)->master_id == master_surface) { this->impactors_information.at(m)->impactor_surfaces->push_back(impactor_surface); break; } } */ ContactRigid::SurfaceToFricCoefMap::iterator it_fc; it_fc = this->friction_coefficient.find(master_surface); if (it_fc != this->friction_coefficient.end()) it_fc->second->addImpactorSurface(impactor_surface); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactRigid::removeMasterSurface(const Surface & master_surface) { AKANTU_DEBUG_IN(); Contact::removeMasterSurface(master_surface); ContactRigid::SurfaceToImpactInfoMap::iterator it_imp; it_imp = this->impactors_information.find(master_surface); AKANTU_DEBUG_ASSERT(it_imp != this->impactors_information.end(), "The master surface: " << master_surface << "couldn't be found and erased in impactors_information map"); delete it_imp->second; this->impactors_information.erase(it_imp); ContactRigid::SurfaceToFricCoefMap::iterator it_fc; it_fc = this->friction_coefficient.find(master_surface); AKANTU_DEBUG_ASSERT(it_fc != this->friction_coefficient.end(), "The master surface: " << master_surface << "couldn't be found and erased in friction_coefficient map"); this->friction_coefficient.erase(it_fc); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactRigid::removeImpactorSurfaceFromMasterSurface(const Surface & impactor_surface, const Surface & master_surface) { AKANTU_DEBUG_IN(); // find in map the impactor information for the given master surface ContactRigid::SurfaceToImpactInfoMap::iterator it_imp; it_imp = this->impactors_information.find(master_surface); AKANTU_DEBUG_ASSERT(it_imp != this->impactors_information.end(), "The master surface: " << master_surface << "couldn't be found and erased in impactors_information map"); std::vector * imp_surfaces = it_imp->second->impactor_surfaces; // find and delete the impactor surface std::vector::iterator it_surf; it_surf = find(imp_surfaces->begin(), imp_surfaces->end(), impactor_surface); AKANTU_DEBUG_ASSERT(it_surf != imp_surfaces->end(), "Couldn't find impactor surface " << impactor_surface << " for the master surface " << master_surface << " and couldn't erase it"); imp_surfaces->erase(it_surf); /* for (UInt m=0; m < this->impactors_information.size(); ++m) { ImpactorInformationPerMaster * impactor_info = this->impactors_information.at(m); if (impactor_info->master_id == master_surface) { for (UInt i=0; i < impactor_info->impactor_surfaces->size(); ++i) { Surface imp_surface = impactor_info->impactor_surfaces->at(i); if (imp_surface == impactor_surface) { impactor_info->impactor_surfaces->erase(impactor_info->impactor_surfaces->begin()+i); break; } } } } */ ContactRigid::SurfaceToFricCoefMap::iterator it_fc; it_fc = this->friction_coefficient.find(master_surface); if (it_fc != this->friction_coefficient.end()) it_fc->second->removeImpactorSurface(impactor_surface); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactRigid::solveContact() { AKANTU_DEBUG_IN(); for(UInt m=0; m < master_surfaces.size(); ++m) { Surface master = this->master_surfaces.at(m); std::stringstream sstr; sstr << id << ":penetration_list"; PenetrationList * penet_list = new PenetrationList(sstr.str()); contact_search->findPenetration(master, *penet_list); solvePenetrationClosestProjection(master, *penet_list); delete penet_list; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /*void ContactRigid::solvePenetration(const PenetrationList & penet_list) { AKANTU_DEBUG_IN(); const UInt dim = ; const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; /// find existing surface element types UInt nb_types = type_list.size(); UInt nb_facet_types = 0; ElementType facet_type[_max_element_type]; for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; if(mesh.getSpatialDimension(type) == spatial_dimension) { ElementType current_facet_type = mesh.getFacetElementType(type); facet_type[nb_facet_types++] = current_facet_type; } } Real * current_position = model.getCurrentPosition().values; Real * displacement = model.getDisplacement().values; Real * increment = model.getIncrement().values; UInt * penetrating_nodes = penet_list.penetrating_nodes.values; for (UInt n=0; n < penet_list.penetrating_nodes.getSize(); ++n) { UInt current_node = penetrating_nodes[n]; // count number of elements penetrated by this node UInt nb_penetrated_facets = 0; ElementType penetrated_type; for (UInt el_type = 0; el_type < nb_facet_types; ++el_type) { ElementType type = facet_type[el_type]; UInt offset_min = penet_list.penetrated_facets_offset[type].get(n); UInt offset_max = penet_list.penetrated_facets_offset[type].get(n+1); nb_penetrated_facets += offset_max - offset_min; penetrated_type = type; } // easy case: node penetrated only one facet if(nb_penetrated_facets == 1) { Real * facets_normals = penet_list.facets_normals[penetrated_type].values; Real * gaps = penet_list.gaps[penetrated_type].values; Real * projected_positions = penet_list.projected_positions[penetrated_type].values; UInt offset_min = penet_list.penetrated_facets_offset[penetrated_type].get(n); for(UInt i=0; i < dim; ++i) { current_position[current_node*dim + i] = projected_positions[offset_min*dim + i]; Real displacement_correction = gaps[offset_min*dim + i] * normals[offset_min*dim + i]; displacement[current_node*dim + i] = displacement[current_node*dim + i] - displacement_correction; increment [current_node*dim + i] = increment [current_node*dim + i] - displacement_correction; } } // if more penetrated facets, find the one from which the impactor node is coming else { } } }*/ /* -------------------------------------------------------------------------- */ void ContactRigid::solvePenetrationClosestProjection(const Surface master, const PenetrationList & penet_list) { AKANTU_DEBUG_IN(); const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; /// find existing surface element types UInt nb_facet_types = 0; ElementType facet_type[_max_element_type]; for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; if(mesh.getSpatialDimension(type) == spatial_dimension) { ElementType current_facet_type = mesh.getFacetElementType(type); facet_type[nb_facet_types++] = current_facet_type; } } for (UInt n=0; n < penet_list.penetrating_nodes.getSize(); ++n) { // find facet for which the gap is the smallest Real min_gap = std::numeric_limits::max(); ElementType penetrated_type; UInt penetrated_facet_offset; for (UInt el_type = 0; el_type < nb_facet_types; ++el_type) { ElementType type = facet_type[el_type]; Real * gaps = penet_list.gaps(type, _not_ghost).storage(); UInt offset_min = penet_list.penetrated_facets_offset(type, _not_ghost)(n); UInt offset_max = penet_list.penetrated_facets_offset(type, _not_ghost)(n+1); for (UInt f = offset_min; f < offset_max; ++f) { if(gaps[f] < min_gap) { min_gap = gaps[f]; penetrated_type = type; penetrated_facet_offset = f; } } } bool is_active = isAlreadyActiveImpactor(master, penet_list, n, penetrated_type, penetrated_facet_offset); if (!is_active) { // correct the position of the impactor projectImpactor(penet_list, n, penetrated_type, penetrated_facet_offset); lockImpactorNode(master, penet_list, n, penetrated_type, penetrated_facet_offset); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ bool ContactRigid::isAlreadyActiveImpactor(const Surface master, const PenetrationList & penet_list, const UInt impactor_index, const ElementType facet_type, const UInt facet_offset) { AKANTU_DEBUG_IN(); bool is_active = false; UInt * penetrating_nodes = penet_list.penetrating_nodes.values; UInt impactor_node = penetrating_nodes[impactor_index]; // find facet normal Real * facets_normals = penet_list.facets_normals(facet_type, _not_ghost).storage(); Real * facet_normal = &facets_normals[facet_offset*spatial_dimension]; Int normal[this->spatial_dimension]; for(UInt i = 0; i < this->spatial_dimension; ++i) normal[i] = Int(floor(facet_normal[i] + 0.5)); // check if this is already in the active impactor node list ContactRigid::SurfaceToImpactInfoMap::iterator it_imp; it_imp = this->impactors_information.find(master); AKANTU_DEBUG_ASSERT(it_imp != this->impactors_information.end(), "The master surface: " << master << "couldn't be found in impactors_information map"); ImpactorInformationPerMaster * impactor_info = it_imp->second; Vector * active_nodes = impactor_info->active_impactor_nodes; Vector * master_normals = impactor_info->master_normals; for (UInt i=0; igetSize(); ++i) { if((*active_nodes)(i) == impactor_node) { UInt count = 0; for (UInt d=0; dspatial_dimension; ++d) { if(Math::are_float_equal((*master_normals)(i,d),normal[d])) count++; } if(count == this->spatial_dimension) is_active = true; } } AKANTU_DEBUG_OUT(); return is_active; } /* -------------------------------------------------------------------------- */ void ContactRigid::projectImpactor(const PenetrationList & penet_list, const UInt impactor_index, const ElementType facet_type, const UInt facet_offset) { AKANTU_DEBUG_IN(); const bool increment_flag = model.getIncrementFlag(); UInt * penetrating_nodes = penet_list.penetrating_nodes.values; Real * facets_normals = penet_list.facets_normals(facet_type, _not_ghost).storage(); Real * gaps = penet_list.gaps(facet_type, _not_ghost).storage(); Real * projected_positions = penet_list.projected_positions(facet_type, _not_ghost).storage(); Real * current_position = model.getCurrentPosition().values; Real * displacement = model.getDisplacement().values; Real * increment = NULL; if(increment_flag) increment = model.getIncrement().values; UInt impactor_node = penetrating_nodes[impactor_index]; for(UInt i=0; i < this->spatial_dimension; ++i) { current_position[impactor_node*spatial_dimension + i] = projected_positions[facet_offset*spatial_dimension + i]; Real displacement_correction = gaps[facet_offset] * facets_normals[facet_offset*spatial_dimension + i]; displacement[impactor_node*spatial_dimension + i] = displacement[impactor_node*spatial_dimension + i] - displacement_correction; if(increment_flag) increment [impactor_node*spatial_dimension + i] = increment [impactor_node*spatial_dimension + i] - displacement_correction; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactRigid::lockImpactorNode(const Surface master, const PenetrationList & penet_list, const UInt impactor_index, const ElementType facet_type, const UInt facet_offset) { AKANTU_DEBUG_IN(); bool init_state_stick = true; // node is in sticking when it gets in contact UInt * penetrating_nodes = penet_list.penetrating_nodes.values; UInt impactor_node = penetrating_nodes[impactor_index]; Real * facets_normals = penet_list.facets_normals(facet_type, _not_ghost).storage(); Real * facet_normal = &facets_normals[facet_offset*spatial_dimension]; Real normal[this->spatial_dimension]; //Int * normal_val = &normal[0]; bool * bound_val = this->model.getBoundary().values; Real * position_val = this->model.getCurrentPosition().values; Real * veloc_val = this->model.getVelocity().values; Real * accel_val = this->model.getAcceleration().values; Real * residual_val = this->model.getResidual().values; for(UInt i = 0; i < this->spatial_dimension; ++i) normal[i] = floor(facet_normal[i] + 0.5); for(UInt i = 0; i < this->spatial_dimension; ++i) { UInt index = impactor_node * spatial_dimension + i; if(!Math::are_float_equal(normal[i],0)) { bound_val[index] = true; veloc_val[index] = 0.; accel_val[index] = 0.; } // if node is in initial stick its tangential velocity has to be zero if(init_state_stick) { veloc_val[index] = 0.; accel_val[index] = 0.; } } ContactRigid::SurfaceToImpactInfoMap::iterator it_imp; it_imp = this->impactors_information.find(master); AKANTU_DEBUG_ASSERT(it_imp != this->impactors_information.end(), "The master surface: " << master << "couldn't be found in impactors_information map"); ImpactorInformationPerMaster * impactor_info = it_imp->second; impactor_info->active_impactor_nodes->push_back(impactor_node); // impactor_info->master_element_offset->push_back(facet_offset); impactor_info->master_element_type->push_back(facet_type); impactor_info->master_normals->push_back(normal); // initial value for stick state when getting in contact bool init_sticking[2]; if (init_state_stick) { init_sticking[0] = true; init_sticking[1] = true; } else { init_sticking[0] = false; init_sticking[1] = false; } impactor_info->node_is_sticking->push_back(init_sticking); Real init_friction_force[this->spatial_dimension]; for(UInt i=0; ispatial_dimension; ++i) { init_friction_force[i] = 0.; } impactor_info->friction_forces->push_back(init_friction_force); impactor_info->stick_positions->push_back(&(position_val[impactor_node*this->spatial_dimension])); impactor_info->residual_forces->push_back(&(residual_val[impactor_node*this->spatial_dimension])); impactor_info->previous_velocities->push_back(&(veloc_val[impactor_node*this->spatial_dimension])); impactor_info->friction_resistances->push_back(0.); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactRigid::avoidAdhesion() { AKANTU_DEBUG_IN(); Real * residual_val = this->model.getResidual().values; bool * bound_val = this->model.getBoundary().values; for(UInt m=0; m < this->master_surfaces.size(); ++m) { Surface master = this->master_surfaces.at(m); ContactRigid::SurfaceToImpactInfoMap::iterator it_imp; it_imp = this->impactors_information.find(master); AKANTU_DEBUG_ASSERT(it_imp != this->impactors_information.end(), "The master surface: " << master << "couldn't be found in impactors_information map"); ImpactorInformationPerMaster * impactor_info = it_imp->second; for (UInt n=0; n < impactor_info->active_impactor_nodes->getSize(); ++n) { UInt current_node = (*impactor_info->active_impactor_nodes)(n); for (UInt i=0; i < spatial_dimension; ++i) { Int direction = Int((*impactor_info->master_normals)(n,i)); Real force = residual_val[current_node * spatial_dimension + i]; if(force * direction > 0.) { bound_val[current_node * spatial_dimension + i] = false; impactor_info->active_impactor_nodes->erase(n); // impactor_info->master_element_offset->erase(n); impactor_info->master_element_type->erase(impactor_info->master_element_type->begin()+n); impactor_info->master_normals->erase(n); impactor_info->node_is_sticking->erase(n); impactor_info->friction_forces->erase(n); impactor_info->stick_positions->erase(n); impactor_info->residual_forces->erase(n); impactor_info->previous_velocities->erase(n); impactor_info->friction_resistances->erase(n); n--; break; } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactRigid::frictionPredictor() { AKANTU_DEBUG_IN(); const Real null_tolerance = std::numeric_limits::epsilon() * 2.; Real * residual_val = this->model.getResidual().values; Real * acceleration_val = this->model.getAcceleration().values; Real * velocity_val = this->model.getVelocity().values; for(UInt m=0; m < this->master_surfaces.size(); ++m) { Surface master = this->master_surfaces.at(m); ContactRigid::SurfaceToImpactInfoMap::iterator it_imp; it_imp = this->impactors_information.find(master); AKANTU_DEBUG_ASSERT(it_imp != this->impactors_information.end(), "The master surface: " << master << "couldn't be found in impactors_information map"); ImpactorInformationPerMaster * impactor_info = it_imp->second; ContactRigid::SurfaceToFricCoefMap::iterator it_fc; it_fc = this->friction_coefficient.find(master); AKANTU_DEBUG_ASSERT(it_fc != this->friction_coefficient.end(), "The master surface: " << master << "couldn't be found in friction_coefficient map"); FrictionCoefficient * fric_coef = it_fc->second; AKANTU_DEBUG_ASSERT(fric_coef != NULL, "There is no friction coefficient defined for master surface " << master); UInt nb_active_impactor_nodes = impactor_info->active_impactor_nodes->getSize(); // compute the friction coefficient for each active impactor node Vector friction_coefficient_values(nb_active_impactor_nodes, 1); Real * friction_coefficient_values_p = friction_coefficient_values.values; fric_coef->computeFrictionCoefficient(friction_coefficient_values); UInt * active_impactor_nodes_val = impactor_info->active_impactor_nodes->storage(); Real * direction_val = impactor_info->master_normals->storage(); bool * node_is_sticking_val = impactor_info->node_is_sticking->storage(); Real * friction_forces_val = impactor_info->friction_forces->storage(); Real * residual_forces_val = impactor_info->residual_forces->storage(); Real * previous_velocities_val = impactor_info->previous_velocities->storage(); Real * friction_resistances_val = impactor_info->friction_resistances->storage(); Vector & nodal_velocity = model.getVelocity(); for (UInt n=0; n < nb_active_impactor_nodes; ++n) { UInt current_node = active_impactor_nodes_val[n]; // save original residual for(UInt i=0; i < spatial_dimension; ++i) { residual_forces_val[n*spatial_dimension+i] = residual_val[current_node*spatial_dimension+i]; } // find friction force mu * normal force Real normal_contact_force = 0.; Real friction_force = 0.; for (UInt i=0; i < spatial_dimension; ++i) { if(!Math::are_float_equal(direction_val[n * this->spatial_dimension + i],0.)) { normal_contact_force = std::abs(residual_val[current_node * this->spatial_dimension + i]); friction_force = friction_coefficient_values_p[n] * normal_contact_force; } } // prakash clifton regularization if (this->prakash) { // compute (|v| + v*)/L Real v_L_term = 0.; for (UInt i=0; i < spatial_dimension; ++i) { if(Math::are_float_equal(direction_val[n * this->spatial_dimension + i],0.)) v_L_term += nodal_velocity(current_node,i) * nodal_velocity(current_node,i); } v_L_term = sqrt(v_L_term); v_L_term += this->ref_velocity; v_L_term /= this->characterstic_length; Real delta_t = this->model.getTimeStep(); // compute friction resistance friction_force *= v_L_term; friction_force += friction_resistances_val[n]/delta_t; friction_force /= (1/delta_t + v_L_term); } + else if (this->dedontney) { + // compute 1/t* + Real v_L_term = 1/ this->t_star; + + Real delta_t = this->model.getTimeStep(); + + // compute friction resistance + friction_force *= v_L_term; + friction_force += friction_resistances_val[n]/delta_t; + friction_force /= (1/delta_t + v_L_term); + } friction_resistances_val[n] = friction_force; // if node is sticking, friction force acts against the residual if (node_is_sticking_val[n*2]) { // compute length of projected residual Real projected_residual = 0.; for (UInt i=0; i < this->spatial_dimension; ++i) { if(Math::are_float_equal(direction_val[n * this->spatial_dimension + i],0.)) { projected_residual += residual_val[current_node * this->spatial_dimension + i] * residual_val[current_node * this->spatial_dimension + i]; } } projected_residual = sqrt(projected_residual); // friction is weaker than residual -> start sliding if (friction_force < projected_residual) { node_is_sticking_val[n*2] = false; node_is_sticking_val[n*2+1] = false; // compute friction vector Real friction_direction[3]; for (UInt i=0; i < this->spatial_dimension; ++i) { if(Math::are_float_equal(direction_val[n * this->spatial_dimension + i],0.)) { friction_direction[i] = residual_val[current_node * this->spatial_dimension + i]; friction_direction[i] /= projected_residual; } else friction_direction[i] = 0.; friction_direction[i] *= -1.; } // add friction force to residual for (UInt i=0; i < this->spatial_dimension; ++i) { Real directional_fric_force = friction_force * friction_direction[i]; friction_forces_val[n*this->spatial_dimension + i] = directional_fric_force; residual_val[current_node * this->spatial_dimension + i] += directional_fric_force; } } // friction is stronger than residual -> stay sticking else { // put residual in friction force and set residual to zero for (UInt i=0; i < this->spatial_dimension; ++i) { if(Math::are_float_equal(direction_val[n * this->spatial_dimension + i],0.)) { friction_forces_val[n*spatial_dimension + i] = residual_val[current_node * this->spatial_dimension + i]; residual_val[current_node * this->spatial_dimension + i] = 0.; } else friction_forces_val[n*spatial_dimension + i] = 0.; friction_forces_val[n*spatial_dimension + i] *= -1.; } } } // end if stick // if node is sliding friction force acts against velocity and // in case there is no velocity it acts against the residual else { Real projected_residual = 0.; Real projected_velocity_magnitude = 0.; for (UInt i=0; i < this->spatial_dimension; ++i) { if(Math::are_float_equal(direction_val[n * this->spatial_dimension + i],0.)) { projected_residual += residual_val[current_node * this->spatial_dimension + i] * residual_val[current_node * this->spatial_dimension + i]; projected_velocity_magnitude += previous_velocities_val[n * this->spatial_dimension + i] * previous_velocities_val[n * this->spatial_dimension + i]; } } projected_residual = sqrt(projected_residual); projected_velocity_magnitude = sqrt(projected_velocity_magnitude); // if no projected velocity nor friction force stronger than projected residual, this is stick if ((projected_velocity_magnitude < null_tolerance) && !(projected_residual > friction_force)) { // put residual in friction force and set residual to zero // set projected velocity and acceleration to zero for (UInt i=0; i < this->spatial_dimension; ++i) { if(Math::are_float_equal(direction_val[n * this->spatial_dimension + i],0.)) { friction_forces_val[n*spatial_dimension + i] = residual_val[current_node * this->spatial_dimension + i]; residual_val [current_node * this->spatial_dimension + i] = 0.; velocity_val [current_node * this->spatial_dimension + i] = 0.; acceleration_val[current_node * this->spatial_dimension + i] = 0.; } else friction_forces_val[n*spatial_dimension + i] = 0.; friction_forces_val[n*spatial_dimension + i] *= -1.; } node_is_sticking_val[n*2] = true; node_is_sticking_val[n*2+1] = true; } // node is really sliding else { UInt index = 0; Real * direction_vector; Real direction_length; // if only velocity is zero (because slider just turned direction) if (projected_velocity_magnitude < null_tolerance) { index = current_node; direction_vector = residual_val; direction_length = projected_residual; } // there is velocity else { index = n; direction_vector = previous_velocities_val; direction_length = projected_velocity_magnitude; } // compute the friction direction Real friction_direction[3]; for (UInt i=0; i < this->spatial_dimension; ++i) { if(Math::are_float_equal(direction_val[n * this->spatial_dimension + i],0.)) { friction_direction[i] = direction_vector[index * this->spatial_dimension + i]; friction_direction[i] /= direction_length; } else friction_direction[i] = 0.; friction_direction[i] *= -1.; } // add friction force to residual for (UInt i=0; i < this->spatial_dimension; ++i) { Real directional_fric_force = friction_force * friction_direction[i]; friction_forces_val[n*this->spatial_dimension + i] = directional_fric_force; residual_val[current_node * this->spatial_dimension + i] += directional_fric_force; } } } // end sliding } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactRigid::frictionCorrector() { AKANTU_DEBUG_IN(); const Real tolerance = std::numeric_limits::epsilon() * 100.; Real * residual_val = this->model.getResidual().values; Real * acceleration_val = this->model.getAcceleration().values; Real * velocity_val = this->model.getVelocity().values; Real time_step = this->model.getTimeStep(); for(UInt m=0; m < this->master_surfaces.size(); ++m) { Surface master = this->master_surfaces.at(m); ContactRigid::SurfaceToImpactInfoMap::iterator it_imp; it_imp = this->impactors_information.find(master); AKANTU_DEBUG_ASSERT(it_imp != this->impactors_information.end(), "The master surface: " << master << "couldn't be found in impactors_information map"); ImpactorInformationPerMaster * impactor_info = it_imp->second; UInt nb_active_impactor_nodes = impactor_info->active_impactor_nodes->getSize(); UInt * active_impactor_nodes_val = impactor_info->active_impactor_nodes->storage(); Real * direction_val = impactor_info->master_normals->storage(); bool * node_is_sticking_val = impactor_info->node_is_sticking->storage(); //Real * friction_forces_val = impactor_info->friction_forces->storage(); //Real * residual_forces_val = impactor_info->residual_forces->storage(); Real * previous_velocities_val = impactor_info->previous_velocities->storage(); for (UInt n=0; n < nb_active_impactor_nodes; ++n) { UInt current_node = active_impactor_nodes_val[n]; // check only sliding impactor nodes if (node_is_sticking_val[n*2]) continue; // compute scalar product of previous velocity with current velocity Real dot_prod_velocities = 0.; Real dot_prod_pred_velocities = 0.; Real current_proj_vel_mag = 0.; for (UInt i=0; i < this->spatial_dimension; ++i) { if(Math::are_float_equal(direction_val[n * this->spatial_dimension + i],0.)) { dot_prod_velocities += velocity_val[current_node * this->spatial_dimension + i] * previous_velocities_val[n * this->spatial_dimension + i]; Real pred_vel = velocity_val[current_node * this->spatial_dimension + i] + 0.5 * time_step * acceleration_val[current_node * this->spatial_dimension + i]; dot_prod_pred_velocities += pred_vel * previous_velocities_val[n * this->spatial_dimension + i]; current_proj_vel_mag += velocity_val[current_node * this->spatial_dimension + i] * velocity_val[current_node * this->spatial_dimension + i]; } } current_proj_vel_mag = sqrt(current_proj_vel_mag); // if velocity has changed sign or has become very small we get into stick if ((dot_prod_velocities < 0.) || (dot_prod_pred_velocities < 0.) || (current_proj_vel_mag < tolerance)) { for (UInt i=0; i < this->spatial_dimension; ++i) { if(Math::are_float_equal(direction_val[n * this->spatial_dimension + i],0.)) { //friction_forces_val[n*spatial_dimension + i] = residual_forces_val[n*spatial_dimension + i]; //residual_forces_val[n * this->spatial_dimension + i] =0.; residual_val [current_node * this->spatial_dimension + i] = 0.; velocity_val [current_node * this->spatial_dimension + i] = 0.; acceleration_val[current_node * this->spatial_dimension + i] = 0.; } //else //friction_forces_val[n*spatial_dimension + i] = 0.; //friction_forces_val[n*spatial_dimension + i] *= -1.; } node_is_sticking_val[n*2] = true; node_is_sticking_val[n*2+1] = true; } } // save current velocities for (UInt n=0; n < nb_active_impactor_nodes; ++n) { UInt current_node = active_impactor_nodes_val[n]; for(UInt i=0; i < spatial_dimension; ++i) { previous_velocities_val[n*spatial_dimension+i] = velocity_val[current_node*spatial_dimension+i]; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactRigid::addRegularizedFriction(const Real & regularizer) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(this->spatial_dimension == 2, "addRegularizedFriction is implemented only for 2D"); Real * residual_val = this->model.getResidual().values; Real * position_val = this->model.getCurrentPosition().values; for(UInt m=0; m < this->master_surfaces.size(); ++m) { Surface master = this->master_surfaces.at(m); // find the impactors information for this master surface ContactRigid::SurfaceToImpactInfoMap::iterator it_imp; it_imp = this->impactors_information.find(master); AKANTU_DEBUG_ASSERT(it_imp != this->impactors_information.end(), "The master surface: " << master << "couldn't be found in impactors_information map"); ImpactorInformationPerMaster * impactor_info = it_imp->second; // find the friction coefficient object for this master ContactRigid::SurfaceToFricCoefMap::iterator it_fc; it_fc = this->friction_coefficient.find(master); AKANTU_DEBUG_ASSERT(it_fc != this->friction_coefficient.end(), "The master surface: " << master << "couldn't be found in friction_coefficient map"); FrictionCoefficient * fric_coef = it_fc->second; AKANTU_DEBUG_ASSERT(fric_coef != NULL, "There is no friction coefficient defined for master surface " << master); UInt nb_active_impactor_nodes = impactor_info->active_impactor_nodes->getSize(); // compute the friction coefficient for each active impactor node Vector friction_coefficient_values(nb_active_impactor_nodes, 1); Real * friction_coefficient_values_p = friction_coefficient_values.values; fric_coef->computeFrictionCoefficient(friction_coefficient_values); UInt * active_impactor_nodes_val = impactor_info->active_impactor_nodes->storage(); Real * direction_val = impactor_info->master_normals->storage(); bool * node_is_sticking_val = impactor_info->node_is_sticking->storage(); Real * friction_forces_val = impactor_info->friction_forces->storage(); Real * stick_positions_val = impactor_info->stick_positions->storage(); for (UInt n=0; n < nb_active_impactor_nodes; ++n) { UInt current_node = active_impactor_nodes_val[n]; Real normal_contact_force = 0.; Real friction_force = 0.; // find friction force mu * normal force for (UInt i=0; ispatial_dimension + i],0.)) { normal_contact_force = std::abs(residual_val[current_node * this->spatial_dimension + i]); friction_force = friction_coefficient_values_p[n] * normal_contact_force; } } // compute F_(i+1)trial Real f_trial = Math::distance_2d(&(position_val[current_node * this->spatial_dimension]), &(stick_positions_val[n * this->spatial_dimension])); f_trial *= regularizer; Real delta_s_vector[this->spatial_dimension]; for (UInt i=0; ispatial_dimension + i] - stick_positions_val[n * this->spatial_dimension + i]; delta_s_vector[i] *= -1.; } // if node is on its own stick position no need to compute friction force // this avoids nan on normalize2 if(Math::norm2(delta_s_vector) < Math::getTolerance()) { node_is_sticking_val[n*2] = true; node_is_sticking_val[n*2+1] = true; continue; } Math::normalize2(delta_s_vector); // check if node is still sticking if (f_trial <= friction_force) { friction_force = f_trial; // sticking position stays unchanged node_is_sticking_val[n*2] = true; node_is_sticking_val[n*2+1] = true; } else { node_is_sticking_val[n*2] = false; node_is_sticking_val[n*2+1] = false; Real delta_s = std::abs(f_trial - friction_force) / regularizer; // friction force stays unchanged for (UInt i=0; ispatial_dimension + i] -= delta_s * delta_s_vector[i]; } } // add friction force to residual for (UInt i=0; i < this->spatial_dimension; ++i) { Real directional_fric_force = friction_force * delta_s_vector[i]; friction_forces_val[n*this->spatial_dimension + i] = directional_fric_force; residual_val[current_node * this->spatial_dimension + i] += directional_fric_force; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactRigid::setStickPositionsToCurrentPositions(const Surface master) { AKANTU_DEBUG_IN(); Real * position_val = this->model.getCurrentPosition().values; // find the impactors information for this master surface ContactRigid::SurfaceToImpactInfoMap::iterator it_imp; it_imp = this->impactors_information.find(master); AKANTU_DEBUG_ASSERT(it_imp != this->impactors_information.end(), "The master surface: " << master << "couldn't be found in impactors_information map"); ImpactorInformationPerMaster * impactor_info = it_imp->second; UInt nb_active_impactor_nodes = impactor_info->active_impactor_nodes->getSize(); UInt * active_impactor_nodes_val = impactor_info->active_impactor_nodes->storage(); Real * stick_positions_val = impactor_info->stick_positions->storage(); for (UInt n=0; n < nb_active_impactor_nodes; ++n) { UInt current_node = active_impactor_nodes_val[n]; // find friction force mu * normal force for (UInt i=0; ispatial_dimension + i] = position_val[current_node*this->spatial_dimension + i]; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactRigid::setFrictionCoefficient(FrictionCoefficient * fric_coefficient) { AKANTU_DEBUG_IN(); Surface master = fric_coefficient->getMasterSurface(); ContactRigid::SurfaceToFricCoefMap::iterator it_fc; it_fc = this->friction_coefficient.find(master); AKANTU_DEBUG_ASSERT(it_fc == this->friction_coefficient.end(), "There is already a friction coefficient for master surface: " << master); this->friction_coefficient[master] = fric_coefficient; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactRigid::removeFrictionCoefficient(const Surface master) { AKANTU_DEBUG_IN(); ContactRigid::SurfaceToFricCoefMap::iterator it_fc; it_fc = this->friction_coefficient.find(master); AKANTU_DEBUG_ASSERT(it_fc != this->friction_coefficient.end(), "There is no friction coefficient for master surface: " << master); this->friction_coefficient.erase(it_fc); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactRigid::getRestartInformation(std::map & map_to_fill, Surface master) { AKANTU_DEBUG_IN(); ContactRigid::SurfaceToImpactInfoMap::iterator it_imp; it_imp = this->impactors_information.find(master); AKANTU_DEBUG_ASSERT(it_imp != this->impactors_information.end(), "The master surface: " << master << "couldn't be found in impactors_information map"); ImpactorInformationPerMaster * impactor_info = it_imp->second; UInt nb_active_nodes = impactor_info->active_impactor_nodes->getSize(); // create vectors Vector * activity_of_nodes = new Vector(this->mesh.getNbNodes(), 1, false); Vector * element_type_of_nodes = new Vector(this->mesh.getNbNodes(), 1, _not_defined); Vector * normals_of_nodes = new Vector(this->mesh.getNbNodes(), this->spatial_dimension, 0.); Vector * sticking_of_nodes = new Vector(this->mesh.getNbNodes(), 2, false); Vector * friction_force_of_nodes = new Vector(this->mesh.getNbNodes(), this->spatial_dimension, 0.); Vector * stick_position_of_nodes = new Vector(this->mesh.getNbNodes(), this->spatial_dimension, 0.); Vector * residual_force_of_nodes = new Vector(this->mesh.getNbNodes(), this->spatial_dimension, 0.); Vector * previous_velocity_of_nodes = new Vector(this->mesh.getNbNodes(), this->spatial_dimension, 0.); Vector * friction_resistances_of_nodes = new Vector(this->mesh.getNbNodes(), 1, 0.); UInt * active_nodes = impactor_info->active_impactor_nodes->storage(); ElementType * element_type = &(*impactor_info->master_element_type)[0]; Real * master_normal = impactor_info->master_normals->storage(); bool * node_stick = impactor_info->node_is_sticking->storage(); Real * friction_force = impactor_info->friction_forces->storage(); Real * stick_position = impactor_info->stick_positions->storage(); Real * residual_force = impactor_info->residual_forces->storage(); Real * previous_velocity = impactor_info->previous_velocities->storage(); Real * friction_resistance = impactor_info->friction_resistances->storage(); for (UInt i=0; ispatial_dimension; ++d, ++master_normal, ++friction_force, ++stick_position, ++residual_force, ++previous_velocity) { (*normals_of_nodes)(*active_nodes,d) = *master_normal; (*friction_force_of_nodes)(*active_nodes,d) = *friction_force; (*stick_position_of_nodes)(*active_nodes,d) = *stick_position; (*residual_force_of_nodes)(*active_nodes,d) = *residual_force; (*previous_velocity_of_nodes)(*active_nodes,d) = *previous_velocity; } for (UInt j=0; j<2; ++j, ++node_stick) { (*sticking_of_nodes)(*active_nodes,j) = *node_stick; } } map_to_fill["active_impactor_nodes"] = activity_of_nodes; map_to_fill["master_element_type"] = element_type_of_nodes; map_to_fill["master_normals"] = normals_of_nodes; map_to_fill["node_is_sticking"] = sticking_of_nodes; map_to_fill["friction_forces"] = friction_force_of_nodes; map_to_fill["stick_positions"] = stick_position_of_nodes; map_to_fill["residual_forces"] = residual_force_of_nodes; map_to_fill["previous_velocities"] = previous_velocity_of_nodes; map_to_fill["friction_resistances"] = friction_resistances_of_nodes; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactRigid::setRestartInformation(std::map & restart_map, Surface master) { AKANTU_DEBUG_IN(); ContactRigid::SurfaceToImpactInfoMap::iterator it_imp; it_imp = this->impactors_information.find(master); AKANTU_DEBUG_ASSERT(it_imp != this->impactors_information.end(), "The master surface: " << master << "couldn't be found in impactors_information map"); ImpactorInformationPerMaster * impactor_info = it_imp->second; std::map < std::string, VectorBase* >::iterator it; it = restart_map.find("active_impactor_nodes"); Vector * ai_nodes = (Vector *)(it->second); if(it == restart_map.end()) { std::cout << "could not find map entry for active impactor nodes" << std::endl; } it = restart_map.find("master_element_type"); Vector * et_nodes = (Vector *)(it->second); if(it == restart_map.end()) { std::cout << "could not find map entry master element type" << std::endl; } it = restart_map.find("master_normals"); Vector * mn_nodes = (Vector *)(it->second); if(it == restart_map.end()) { std::cout << "could not find map entry for master normals" << std::endl; } it = restart_map.find("node_is_sticking"); Vector * is_nodes = (Vector *)(it->second); if(it == restart_map.end()) { std::cout << "could not find map entry node is sticking" << std::endl; } it = restart_map.find("friction_forces"); Vector * ff_nodes = (Vector *)(it->second); if(it == restart_map.end()) { std::cout << "could not find map entry friction forces" << std::endl; } it = restart_map.find("stick_positions"); Vector * sp_nodes = (Vector *)(it->second); if(it == restart_map.end()) { std::cout << "could not find map entry stick positions" << std::endl; } it = restart_map.find("residual_forces"); Vector * rf_nodes = (Vector *)(it->second); if(it == restart_map.end()) { std::cout << "could not find map entry residual forces" << std::endl; } it = restart_map.find("previous_velocities"); Vector * pv_nodes = (Vector *)(it->second); if(it == restart_map.end()) { std::cout << "could not find map entry previous velocities" << std::endl; } it = restart_map.find("friction_resistances"); Vector * fr_nodes = (Vector *)(it->second); if(it == restart_map.end()) { std::cout << "could not find map entry friction resistances" << std::endl; } UInt nb_nodes = this->mesh.getNbNodes(); for (UInt i=0; iactive_impactor_nodes->push_back(i); impactor_info->master_element_type->push_back((*et_nodes)(i)); // get master_normals and friction_forces information Real normal[this->spatial_dimension]; Real friction[this->spatial_dimension]; Real position[this->spatial_dimension]; Real residual[this->spatial_dimension]; Real velocity[this->spatial_dimension]; for (UInt d=0; dspatial_dimension; ++d) { normal[d] = (*mn_nodes)(i,d); friction[d] = (*ff_nodes)(i,d); position[d] = (*sp_nodes)(i,d); residual[d] = (*rf_nodes)(i,d); velocity[d] = (*pv_nodes)(i,d); } impactor_info->master_normals->push_back(normal); impactor_info->friction_forces->push_back(friction); impactor_info->stick_positions->push_back(position); impactor_info->residual_forces->push_back(residual); impactor_info->previous_velocities->push_back(velocity); impactor_info->friction_resistances->push_back((*fr_nodes)(i)); // get node_is_sticking information bool stick[2]; for (UInt j=0; j<2; ++j) { stick[j] = (*is_nodes)(i,j); } impactor_info->node_is_sticking->push_back(stick); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void ContactRigid::setSimplifiedPrakashCliftonFriction(Real v_star, Real length) { this->prakash = true; this->ref_velocity = v_star; this->characterstic_length = length; } /* -------------------------------------------------------------------------- */ void ContactRigid::unsetSimplifiedPrakashCliftonFriction() { this->prakash = false; + this->dedontney = false; +} + +/* -------------------------------------------------------------------------- */ +void ContactRigid::setSimplifiedPrakashCliftonFriction(Real t_star) { + this->dedontney = true; + this->t_star = t_star; } /* -------------------------------------------------------------------------- */ void ContactRigid::setPrakashCliftonToSteadyState(const Surface master) { AKANTU_DEBUG_IN(); Real * residual_val = this->model.getResidual().values; // find the impactors information for this master surface ContactRigid::SurfaceToImpactInfoMap::iterator it_imp; it_imp = this->impactors_information.find(master); AKANTU_DEBUG_ASSERT(it_imp != this->impactors_information.end(), "The master surface: " << master << "couldn't be found in impactors_information map"); ImpactorInformationPerMaster * impactor_info = it_imp->second; ContactRigid::SurfaceToFricCoefMap::iterator it_fc; it_fc = this->friction_coefficient.find(master); AKANTU_DEBUG_ASSERT(it_fc != this->friction_coefficient.end(), "The master surface: " << master << "couldn't be found in friction_coefficient map"); FrictionCoefficient * fric_coef = it_fc->second; AKANTU_DEBUG_ASSERT(fric_coef != NULL, "There is no friction coefficient defined for master surface " << master); UInt nb_active_impactor_nodes = impactor_info->active_impactor_nodes->getSize(); UInt * active_impactor_nodes_val = impactor_info->active_impactor_nodes->storage(); Real * direction_val = impactor_info->master_normals->storage(); Real * friction_resistances_val = impactor_info->friction_resistances->storage(); // compute the friction coefficient for each active impactor node Vector friction_coefficient_values(nb_active_impactor_nodes, 1); Real * friction_coefficient_values_p = friction_coefficient_values.values; fric_coef->computeFrictionCoefficient(friction_coefficient_values); for (UInt n=0; n < nb_active_impactor_nodes; ++n) { UInt current_node = active_impactor_nodes_val[n]; Real normal_contact_force = 0.; Real friction_force = 0.; for (UInt i=0; i < spatial_dimension; ++i) { if(!Math::are_float_equal(direction_val[n * this->spatial_dimension + i],0.)) { normal_contact_force = std::abs(residual_val[current_node * this->spatial_dimension + i]); friction_force = friction_coefficient_values_p[n] * normal_contact_force; } } friction_resistances_val[n] = friction_force; } AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/src/model/solid_mechanics/contact/contact_rigid.hh b/src/model/solid_mechanics/contact/contact_rigid.hh index 188d2c00c..26cded4c3 100644 --- a/src/model/solid_mechanics/contact/contact_rigid.hh +++ b/src/model/solid_mechanics/contact/contact_rigid.hh @@ -1,245 +1,248 @@ /** * @file contact_rigid.hh * * @author David Simon Kammer * * @date Wed Jan 19 14:45:41 2011 * * @brief Structure that solves contact for for a rigid master surface without * friction within an explicit time scheme * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_CONTACT_RIGID_HH__ #define __AKANTU_CONTACT_RIGID_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "contact.hh" #include "friction_coefficient.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ class ContactRigid : public Contact { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: ContactRigid(const SolidMechanicsModel & model, const ContactType & type, const ContactID & id = "contact", const MemoryID & memory_id = 0); virtual ~ContactRigid(); /* ------------------------------------------------------------------------ */ /* Active Contact Nodes Class (per Master Surface) */ /* ------------------------------------------------------------------------ */ public: class ImpactorInformationPerMaster { public: ImpactorInformationPerMaster(const Surface master_id, const UInt spatial_dimension); virtual ~ImpactorInformationPerMaster(); public: /// master surface id Surface master_id; /// impactor surfaces std::vector * impactor_surfaces; /// list of active impactor nodes Vector * active_impactor_nodes; // the offset of the associated master surface elemet //Vector * master_element_offset; /// the element type of the associated master surface element std::vector * master_element_type; /// the normal to the master surface element for each active impactor node Vector * master_normals; /// show if node is sticking Vector * node_is_sticking; /// friction force for each Vector * friction_forces; /// stick position for regularized friction Vector * stick_positions; /// residual forces without friction force Vector * residual_forces; /// velocities before predictor computation Vector * previous_velocities; /// friction resistance Vector * friction_resistances; }; typedef std::map SurfaceToImpactInfoMap; typedef std::map SurfaceToFricCoefMap; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// solve the contact virtual void solveContact(); /// avoid adhesion by delocking contact nodes that have tensile contact force virtual void avoidAdhesion(); /// alternative way for addFriction and addSticking virtual void addRegularizedFriction(const Real & regularizer); /// friction predictor virtual void frictionPredictor(); /// friction corrector virtual void frictionCorrector(); /// reset stick positions to current positions virtual void setStickPositionsToCurrentPositions(const Surface master); /// add a new master surface virtual void addMasterSurface(const Surface & master_surface); /// add an impactor surface to a master surface virtual void addImpactorSurfaceToMasterSurface(const Surface & impactor_surface, const Surface & master_surface); /// remove a master surface virtual void removeMasterSurface(const Surface & master_surface); /// remove an impactor surface from master surface virtual void removeImpactorSurfaceFromMasterSurface(const Surface & impactor_surface, const Surface & master_surface); /// set a friction coefficient for a given master surface virtual void setFrictionCoefficient(FrictionCoefficient * fric_coefficient); /// remove a friction coefficient for a given master surface virtual void removeFrictionCoefficient(const Surface master); /// put contact information into a map which can be used for restart virtual void getRestartInformation(std::map & map_to_fill, Surface master); /// put contact information into a map which can be used for restart virtual void setRestartInformation(std::map & restart_map, Surface master); /// compute friction based on the simplified prakash clifton regularization (cochard & rice 2000) void setSimplifiedPrakashCliftonFriction(Real v_star, Real length); + void setSimplifiedPrakashCliftonFriction(Real t_star); void unsetSimplifiedPrakashCliftonFriction(); void setPrakashCliftonToSteadyState(const Surface master); /// function to print the contain of the class //virtual void printself(std::ostream & stream, int indent = 0) const; private: /// solve penetration of penetrating nodes //void solvePenetration(const PenetrationList & penet_list); /// solve penetration to the closest facet void solvePenetrationClosestProjection(const Surface master, const PenetrationList & penet_list); /// find if the penetrated node is not already in list for this master with this normal bool isAlreadyActiveImpactor(const Surface master, const PenetrationList & penet_list, const UInt impactor_index, const ElementType facet_type, const UInt facet_offset); /// projects the impactor to the projected position void projectImpactor(const PenetrationList & penet_list, const UInt impactor_index, const ElementType facet_type, const UInt facet_offset); void lockImpactorNode(const Surface master, const PenetrationList & penet_list, const UInt impactor_index, const ElementType facet_type, const UInt facet_offset); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the vector full of impactor information objects AKANTU_GET_MACRO(ImpactorsInformation, impactors_information, const SurfaceToImpactInfoMap &); /// get spatial dimension AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// spatial dimension of mesh UInt spatial_dimension; /// the mesh const Mesh & mesh; // simplified prakash clifton regularization bool prakash; + bool dedontney; // parameters needed for fomulation of cochard and rice 2000 Real ref_velocity; Real characterstic_length; + Real t_star; /// list of impactor nodes info for each master surface SurfaceToImpactInfoMap impactors_information; /// list of friction coefficient objects for each master surface SurfaceToFricCoefMap friction_coefficient; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "contact_rigid_inline_impl.cc" /// standard output stream operator //inline std::ostream & operator <<(std::ostream & stream, const ContactRigid & _this) //{ // _this.printself(stream); // return stream; //} __END_AKANTU__ #endif /*__AKANTU_CONTACT_RIGID_HH__ */ diff --git a/src/model/solid_mechanics/material.cc b/src/model/solid_mechanics/material.cc index b6c44639a..96ae60c41 100644 --- a/src/model/solid_mechanics/material.cc +++ b/src/model/solid_mechanics/material.cc @@ -1,988 +1,987 @@ /** * @file material.cc * * @author Marco Vocialta * @author Nicolas Richart * * @date Tue Jul 27 18:15:37 2010 * * @brief Implementation of the common part of the material class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "solid_mechanics_model.hh" #include "sparse_matrix.hh" #include "dof_synchronizer.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ Material::Material(SolidMechanicsModel & model, const ID & id) : Memory(model.getMemoryID()), id(id), name(""), model(&model), stress("stress", id), strain("strain", id), element_filter("element_filter", id), // potential_energy_vector(false), potential_energy("potential_energy", id), is_non_local(false), interpolation_inverse_coordinates("interpolation inverse coordinates", id), interpolation_points_matrices("interpolation points matrices", id), is_init(false) { AKANTU_DEBUG_IN(); registerParam("rho", rho, 0., ParamAccessType(_pat_parsable | _pat_modifiable), "Density"); registerParam("id", this->id, _pat_readable); registerParam("name", name, std::string(), ParamAccessType(_pat_parsable | _pat_readable)); spatial_dimension = this->model->getSpatialDimension(); /// allocate strain stress for local elements initInternalVector(strain, spatial_dimension * spatial_dimension); initInternalVector(stress, spatial_dimension * spatial_dimension); /// for each connectivity types allocate the element filer array of the material initInternalVector(element_filter, 1, true); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Material::~Material() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ bool Material::parseParam(const std::string & key, const std::string & value, __attribute__ ((unused)) const ID & id) { try { params.parseParam(key, value); } catch(...) { return false; } return true; } /* -------------------------------------------------------------------------- */ void Material::initMaterial() { AKANTU_DEBUG_IN(); resizeInternalVector(stress); resizeInternalVector(strain); is_init = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::initInternalVector(ByElementTypeVector & vect, UInt nb_component, bool temporary, ElementKind element_kind) { AKANTU_DEBUG_IN(); model->getFEM().getMesh().initByElementTypeVector(vect, nb_component, spatial_dimension, false, element_kind); if(!temporary) registerInternal(vect); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template<> void Material::registerInternal(ByElementTypeVector & vect) { internal_vectors_real[vect.getID()] = &vect; } template<> void Material::registerInternal(ByElementTypeVector & vect) { internal_vectors_uint[vect.getID()] = &vect; } /* -------------------------------------------------------------------------- */ template void Material::resizeInternalVector(ByElementTypeVector & by_el_type_vect, ElementKind element_kind) const { AKANTU_DEBUG_IN(); FEM * fem = & model->getFEM(); if (element_kind == _ek_cohesive) fem = & model->getFEM("CohesiveFEM"); const Mesh & mesh = fem->getMesh(); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, element_kind); Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, element_kind); for(; it != end; ++it) { const Vector & elem_filter = element_filter(*it, gt); UInt nb_element = elem_filter.getSize(); UInt nb_quadrature_points = fem->getNbQuadraturePoints(*it, gt); UInt new_size = nb_element * nb_quadrature_points; Vector & vect = by_el_type_vect(*it, gt); vect.resize(new_size); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the residual by assembling @f$\int_{e} \sigma_e \frac{\partial * \varphi}{\partial X} dX @f$ * * @param[in] displacements nodes displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::updateResidual(GhostType ghost_type) { AKANTU_DEBUG_IN(); computeAllStresses(ghost_type); assembleResidual(ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::assembleResidual(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Vector & residual = const_cast &>(model->getResidual()); Mesh & mesh = model->getFEM().getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { const Vector & shapes_derivatives = model->getFEM().getShapesDerivatives(*it, ghost_type); Vector & elem_filter = element_filter(*it, ghost_type); UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(*it, ghost_type); UInt nb_element = elem_filter.getSize(); /// compute @f$\sigma \frac{\partial \varphi}{\partial X}@f$ by @f$\mathbf{B}^t \mathbf{\sigma}_q@f$ Vector * sigma_dphi_dx = new Vector(nb_element*nb_quadrature_points, size_of_shapes_derivatives, "sigma_x_dphi_/_dX"); Real * shapesd = shapes_derivatives.storage(); Real * shapesd_val; UInt * elem_filter_val = elem_filter.storage(); Vector * shapesd_filtered = new Vector(nb_element*nb_quadrature_points, size_of_shapes_derivatives, "filtered shapesd"); Real * shapesd_filtered_val = shapesd_filtered->values; for (UInt el = 0; el < nb_element; ++el) { shapesd_val = shapesd + elem_filter_val[el] * size_of_shapes_derivatives * nb_quadrature_points; memcpy(shapesd_filtered_val, shapesd_val, size_of_shapes_derivatives * nb_quadrature_points * sizeof(Real)); shapesd_filtered_val += size_of_shapes_derivatives * nb_quadrature_points; } Vector & stress_vect = stress(*it, ghost_type); // Vector::iterator sigma = stress_vect.begin(spatial_dimension, spatial_dimension); // Vector::iterator sigma_end = stress_vect.end(spatial_dimension, spatial_dimension); // Vector::iterator nabla_B = shapesd_filtered->begin(nb_nodes_per_element, spatial_dimension); // Vector::iterator sigma_dphi_dx_it = sigma_dphi_dx->begin(nb_nodes_per_element, spatial_dimension); // for (; sigma != sigma_end; ++sigma, ++nabla_B, ++sigma_dphi_dx_it) { // sigma_dphi_dx_it->mul(*nabla_B, *sigma); // } Math::matrix_matrixt(nb_nodes_per_element, spatial_dimension, spatial_dimension, *shapesd_filtered, stress_vect, *sigma_dphi_dx); delete shapesd_filtered; /** * compute @f$\int \sigma * \frac{\partial \varphi}{\partial X}dX@f$ by @f$ \sum_q \mathbf{B}^t * \mathbf{\sigma}_q \overline w_q J_q@f$ */ Vector * int_sigma_dphi_dx = new Vector(nb_element, nb_nodes_per_element * spatial_dimension, "int_sigma_x_dphi_/_dX"); model->getFEM().integrate(*sigma_dphi_dx, *int_sigma_dphi_dx, size_of_shapes_derivatives, *it, ghost_type, &elem_filter); delete sigma_dphi_dx; /// assemble model->getFEM().assembleVector(*int_sigma_dphi_dx, residual, model->getDOFSynchronizer().getLocalDOFEquationNumbers(), residual.getNbComponent(), *it, ghost_type, &elem_filter, -1); delete int_sigma_dphi_dx; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the stress from the strain * * @param[in] current_position nodes postition + displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::computeAllStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); resizeInternalVector(stress); resizeInternalVector(strain); UInt spatial_dimension = model->getSpatialDimension(); Mesh::type_iterator it = model->getFEM().getMesh().firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = model->getFEM().getMesh().lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { Vector & elem_filter = element_filter(*it, ghost_type); Vector & strain_vect = strain(*it, ghost_type); /// compute @f$\nabla u@f$ model->getFEM().gradientOnQuadraturePoints(model->getDisplacement(), strain_vect, spatial_dimension, *it, ghost_type, &elem_filter); /// compute @f$\mathbf{\sigma}_q@f$ from @f$\nabla u@f$ computeStress(*it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::setToSteadyState(GhostType ghost_type) { AKANTU_DEBUG_IN(); const Vector & displacement = model->getDisplacement(); resizeInternalVector(strain); UInt spatial_dimension = model->getSpatialDimension(); Mesh::type_iterator it = model->getFEM().getMesh().firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = model->getFEM().getMesh().lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { Vector & elem_filter = element_filter(*it, ghost_type); Vector & strain_vect = strain(*it, ghost_type); /// compute @f$\nabla u@f$ model->getFEM().gradientOnQuadraturePoints(displacement, strain_vect, spatial_dimension, *it, ghost_type, &elem_filter); setToSteadyState(*it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the stiffness matrix by assembling @f$\int_{\omega} B^t \times D * \times B d\omega @f$ * * @param[in] current_position nodes postition + displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::assembleStiffnessMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Mesh & mesh = model->getFEM().getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { switch(spatial_dimension) { case 1: { assembleStiffnessMatrix<1>(*it, ghost_type); break; } case 2: { assembleStiffnessMatrix<2>(*it, ghost_type); break; } case 3: { assembleStiffnessMatrix<3>(*it, ghost_type); break; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); SparseMatrix & K = const_cast(model->getStiffnessMatrix()); const Vector & shapes_derivatives = model->getFEM().getShapesDerivatives(type,ghost_type); Vector & elem_filter = element_filter(type, ghost_type); Vector & strain_vect = strain(type, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(type, ghost_type); strain_vect.resize(nb_quadrature_points * nb_element); model->getFEM().gradientOnQuadraturePoints(model->getDisplacement(), strain_vect, dim, type, ghost_type, &elem_filter); UInt tangent_size = getTangentStiffnessVoigtSize(dim); Vector * tangent_stiffness_matrix = new Vector(nb_element*nb_quadrature_points, tangent_size * tangent_size, "tangent_stiffness_matrix"); tangent_stiffness_matrix->clear(); computeTangentModuli(type, *tangent_stiffness_matrix, ghost_type); Vector * shapes_derivatives_filtered = new Vector(nb_element * nb_quadrature_points, dim * nb_nodes_per_element, "shapes derivatives filtered"); Vector::const_iterator shapes_derivatives_it = shapes_derivatives.begin(spatial_dimension, nb_nodes_per_element); Vector::iterator shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(spatial_dimension, nb_nodes_per_element); UInt * elem_filter_val = elem_filter.storage(); for (UInt e = 0; e < nb_element; ++e, ++elem_filter_val) for (UInt q = 0; q < nb_quadrature_points; ++q, ++shapes_derivatives_filtered_it) *shapes_derivatives_filtered_it = shapes_derivatives_it[*elem_filter_val * nb_quadrature_points + q]; /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_d_b_size = dim * nb_nodes_per_element; Vector * bt_d_b = new Vector(nb_element * nb_quadrature_points, bt_d_b_size * bt_d_b_size, "B^t*D*B"); types::RMatrix B(tangent_size, dim * nb_nodes_per_element); types::RMatrix Bt_D(dim * nb_nodes_per_element, tangent_size); shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(nb_nodes_per_element, spatial_dimension); Vector::iterator Bt_D_B_it = bt_d_b->begin(dim*nb_nodes_per_element, dim*nb_nodes_per_element); Vector::iterator D_it = tangent_stiffness_matrix->begin(tangent_size, tangent_size); Vector::iterator D_end = tangent_stiffness_matrix->end (tangent_size, tangent_size); for(; D_it != D_end; ++D_it, ++Bt_D_B_it, ++shapes_derivatives_filtered_it) { types::RMatrix & D = *D_it; types::RMatrix & Bt_D_B = *Bt_D_B_it; transferBMatrixToSymVoigtBMatrix(*shapes_derivatives_filtered_it, B, nb_nodes_per_element); Bt_D.mul(B, D); Bt_D_B.mul(Bt_D, B); } delete tangent_stiffness_matrix; delete shapes_derivatives_filtered; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Vector * K_e = new Vector(nb_element, bt_d_b_size * bt_d_b_size, "K_e"); model->getFEM().integrate(*bt_d_b, *K_e, bt_d_b_size * bt_d_b_size, type, ghost_type, &elem_filter); delete bt_d_b; model->getFEM().assembleMatrix(*K_e, K, spatial_dimension, type, ghost_type, &elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computeAllStressesFromTangentModuli(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Mesh & mesh = model->getFEM().getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { switch(spatial_dimension) { case 1: { computeAllStressesFromTangentModuli<1>(*it, ghost_type); break; } case 2: { computeAllStressesFromTangentModuli<2>(*it, ghost_type); break; } case 3: { computeAllStressesFromTangentModuli<3>(*it, ghost_type); break; } } } AKANTU_DEBUG_OUT(); } - /* -------------------------------------------------------------------------- */ template void Material::computeAllStressesFromTangentModuli(const ElementType & type, - GhostType ghost_type) { + GhostType ghost_type) { AKANTU_DEBUG_IN(); const Vector & shapes_derivatives = model->getFEM().getShapesDerivatives(type, ghost_type); Vector & elem_filter = element_filter(type, ghost_type); Vector & strain_vect = strain(type, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(type, ghost_type); const Vector & connectivity = model->getFEM().getMesh().getConnectivity(type, ghost_type); strain_vect.resize(nb_quadrature_points * nb_element); Vector & u = model->getDisplacement(); model->getFEM().gradientOnQuadraturePoints(u, strain_vect, dim, type, ghost_type, &elem_filter); UInt tangent_moduli_size = getTangentStiffnessVoigtSize(dim); Vector * tangent_moduli_tensors = new Vector(nb_element*nb_quadrature_points, tangent_moduli_size * tangent_moduli_size, "tangent_moduli_tensors"); tangent_moduli_tensors->clear(); computeTangentModuli(type, *tangent_moduli_tensors, ghost_type); Vector * shapes_derivatives_filtered = new Vector(nb_element * nb_quadrature_points, dim * nb_nodes_per_element, "shapes derivatives filtered"); Vector::const_iterator shapes_derivatives_it = shapes_derivatives.begin(spatial_dimension, nb_nodes_per_element); Vector::iterator shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(spatial_dimension, nb_nodes_per_element); Vector filtered_u(nb_element, nb_nodes_per_element * spatial_dimension); UInt * elem_filter_val = elem_filter.storage(); Real * filtered_u_it = filtered_u.storage(); for (UInt e = 0; e < nb_element; ++e, ++elem_filter_val) { UInt el = *elem_filter_val; for (UInt n = 0; n < nb_nodes_per_element; ++n) { for (UInt s = 0; s < spatial_dimension; ++s, ++filtered_u_it) { - *filtered_u_it = u(connectivity(el,n), s); + *filtered_u_it = u(connectivity(el,n), s); } } for (UInt q = 0; q < nb_quadrature_points; ++q, ++shapes_derivatives_filtered_it) *shapes_derivatives_filtered_it = shapes_derivatives_it[el * nb_quadrature_points + q]; } /// compute @f$\mathbf{D} \mathbf{B} \mathbf{u}@f$ shapes_derivatives_filtered_it = shapes_derivatives_filtered->begin(nb_nodes_per_element, spatial_dimension); Vector::iterator D_it = tangent_moduli_tensors->begin(tangent_moduli_size, - tangent_moduli_size); + tangent_moduli_size); Vector::iterator sigma_it = stress(type, ghost_type).begin(spatial_dimension, - spatial_dimension); + spatial_dimension); Vector::iterator u_it = filtered_u.begin(spatial_dimension * nb_nodes_per_element); types::RMatrix B(tangent_moduli_size, spatial_dimension * nb_nodes_per_element); types::RVector Bu(tangent_moduli_size); types::RVector DBu(tangent_moduli_size); for (UInt e = 0; e < nb_element; ++e, ++u_it) { for (UInt q = 0; q < nb_quadrature_points; ++q, ++D_it, ++shapes_derivatives_filtered_it, ++sigma_it) { types::RVector & u = *u_it; types::RMatrix & sigma = *sigma_it; types::RMatrix & D = *D_it; transferBMatrixToSymVoigtBMatrix(*shapes_derivatives_filtered_it, B, nb_nodes_per_element); Bu.mul(B, u); DBu.mul(D, Bu); // Voigt notation to full symmetric tensor for (UInt i = 0; i < dim; ++i) sigma(i, i) = DBu(i); if(dim == 2) { sigma(0,1) = sigma(1,0) = DBu(2); } else if(dim == 3) { sigma(1,2) = sigma(2,1) = DBu(3); sigma(0,2) = sigma(2,0) = DBu(4); sigma(0,1) = sigma(1,0) = DBu(5); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergy(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); if(ghost_type != _not_ghost) return; Real * epot = potential_energy(el_type, ghost_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); computePotentialEnergyOnQuad(grad_u, sigma, *epot); epot++; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergyByElement() { AKANTU_DEBUG_IN(); Mesh & mesh = model->getFEM().getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension); for(; it != last_type; ++it) { if(!potential_energy.exists(*it, _not_ghost)) { UInt nb_element = element_filter(*it, _not_ghost).getSize(); UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(*it, _not_ghost); potential_energy.alloc(nb_element * nb_quadrature_points, 1, *it, _not_ghost); } computePotentialEnergy(*it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergyByElement(ElementType type, UInt index, types::RVector & epot_on_quad_points){ - Vector::iterator strain_it = - this->strain(type).begin(spatial_dimension, - spatial_dimension); - Vector::iterator strain_end = - this->strain(type).begin(spatial_dimension, - spatial_dimension); - Vector::iterator stress_it = - this->stress(type).begin(spatial_dimension, - spatial_dimension); - + Vector::iterator strain_it = + this->strain(type).begin(spatial_dimension, + spatial_dimension); + Vector::iterator strain_end = + this->strain(type).begin(spatial_dimension, + spatial_dimension); + Vector::iterator stress_it = + this->stress(type).begin(spatial_dimension, + spatial_dimension); + UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(type); strain_it += index*nb_quadrature_points; strain_end += (index+1)*nb_quadrature_points; stress_it += index*nb_quadrature_points; - + Real * epot_quad = epot_on_quad_points.storage(); - for(;strain_it != strain_end; ++strain_it, ++stress_it, ++ epot_quad) { - types::RMatrix & __attribute__((unused)) grad_u = *strain_it; + for(;strain_it != strain_end; ++strain_it, ++stress_it, ++ epot_quad) { + types::RMatrix & __attribute__((unused)) grad_u = *strain_it; types::RMatrix & __attribute__((unused)) sigma = *stress_it; computePotentialEnergyOnQuad(grad_u,sigma,*epot_quad); } } /* -------------------------------------------------------------------------- */ Real Material::getPotentialEnergy() { AKANTU_DEBUG_IN(); Real epot = 0.; computePotentialEnergyByElement(); /// integrate the potential energy for each type of elements Mesh & mesh = model->getFEM().getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension); for(; it != last_type; ++it) { epot += model->getFEM().integrate(potential_energy(*it, _not_ghost), *it, _not_ghost, &element_filter(*it, _not_ghost)); } AKANTU_DEBUG_OUT(); return epot; } /* -------------------------------------------------------------------------- */ Real Material::getPotentialEnergy(ElementType & type, UInt index) { AKANTU_DEBUG_IN(); Real epot = 0.; types::RVector epot_on_quad_points(model->getFEM().getNbQuadraturePoints(type)); - + computePotentialEnergyByElement(type,index,epot_on_quad_points); epot = model->getFEM().integrate(epot_on_quad_points, type, element_filter(type)(index)); AKANTU_DEBUG_OUT(); return epot; } /* -------------------------------------------------------------------------- */ Real Material::getEnergy(std::string type) { AKANTU_DEBUG_IN(); if(type == "potential") return getPotentialEnergy(); AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ Real Material::getEnergy(std::string energy_id, ElementType type, UInt index) { AKANTU_DEBUG_IN(); if(energy_id == "potential") return getPotentialEnergy(type,index); AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ void Material::computeQuadraturePointsCoordinates(ByElementTypeReal & quadrature_points_coordinates, const GhostType & ghost_type) const { AKANTU_DEBUG_IN(); const Mesh & mesh = model->getFEM().getMesh(); Vector nodes_coordinates(mesh.getNodes(), true); nodes_coordinates += model->getDisplacement(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { const Vector & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_tot_quad = model->getFEM().getNbQuadraturePoints(*it, ghost_type) * nb_element; Vector & quads = quadrature_points_coordinates(*it, ghost_type); quads.resize(nb_tot_quad); model->getFEM().interpolateOnQuadraturePoints(nodes_coordinates, quads, spatial_dimension, *it, ghost_type, &elem_filter); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::initElementalFieldInterpolation(ByElementTypeReal & interpolation_points_coordinates) { AKANTU_DEBUG_IN(); const Mesh & mesh = model->getFEM().getMesh(); ByElementTypeReal quadrature_points_coordinates("quadrature_points_coordinates_tmp_nl", id); mesh.initByElementTypeVector(quadrature_points_coordinates, spatial_dimension, 0); computeQuadraturePointsCoordinates(quadrature_points_coordinates, _not_ghost); computeQuadraturePointsCoordinates(quadrature_points_coordinates, _ghost); Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator last = mesh.lastType(spatial_dimension); for (; it != last; ++it) { UInt nb_element = mesh.getNbElement(*it); if (nb_element == 0) continue; ElementType type = *it; #define AKANTU_INIT_INTERPOLATE_ELEMENTAL_FIELD(type) \ initElementalFieldInterpolation(quadrature_points_coordinates(type), \ interpolation_points_coordinates(type)) AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(AKANTU_INIT_INTERPOLATE_ELEMENTAL_FIELD); #undef AKANTU_INIT_INTERPOLATE_ELEMENTAL_FIELD } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::initElementalFieldInterpolation(const Vector & quad_coordinates, const Vector & interpolation_points_coordinates) { AKANTU_DEBUG_IN(); UInt size_inverse_coords = getSizeElementalFieldInterpolationCoodinates(); Vector & elem_fil = element_filter(type); UInt nb_element = elem_fil.getSize(); UInt nb_quad_per_element = model->getFEM().getNbQuadraturePoints(type); UInt nb_interpolation_points = interpolation_points_coordinates.getSize(); AKANTU_DEBUG_ASSERT(nb_interpolation_points % nb_element == 0, "Can't interpolate elemental field on elements, the coordinates vector has a wrong size"); UInt nb_interpolation_points_per_elem = nb_interpolation_points / nb_element; if(!interpolation_inverse_coordinates.exists(type)) interpolation_inverse_coordinates.alloc(nb_element, size_inverse_coords*size_inverse_coords, type); if(!interpolation_points_matrices.exists(type)) interpolation_points_matrices.alloc(nb_element, nb_interpolation_points_per_elem * size_inverse_coords, type); Vector & interp_inv_coord = interpolation_inverse_coordinates(type); Vector & interp_points_mat = interpolation_points_matrices(type); types::RMatrix quad_coord_matrix(size_inverse_coords, size_inverse_coords); Vector::const_iterator quad_coords_it = quad_coordinates.begin_reinterpret(nb_quad_per_element, spatial_dimension, nb_element); Vector::const_iterator points_coords_it = interpolation_points_coordinates.begin_reinterpret(nb_interpolation_points_per_elem, spatial_dimension, nb_element); Vector::iterator inv_quad_coord_it = interp_inv_coord.begin(size_inverse_coords, size_inverse_coords); Vector::iterator inv_points_mat_it = interp_points_mat.begin(nb_interpolation_points_per_elem, size_inverse_coords); /// loop over the elements of the current material and element type for (UInt el = 0; el < nb_element; ++el) { /// matrix containing the quadrature points coordinates const types::RMatrix & quad_coords = *quad_coords_it; /// matrix to store the matrix inversion result types::RMatrix & inv_quad_coord_matrix = *inv_quad_coord_it; /// insert the quad coordinates in a matrix compatible with the interpolation buildElementalFieldInterpolationCoodinates(quad_coords, quad_coord_matrix); /// invert the interpolation matrix inv_quad_coord_matrix.inverse(quad_coord_matrix); /// matrix containing the interpolation points coordinates const types::RMatrix & points_coords = *points_coords_it; /// matrix to store the interpolation points coordinates /// compatible with these functions types::RMatrix & inv_points_coord_matrix = *inv_points_mat_it; /// insert the quad coordinates in a matrix compatible with the interpolation buildElementalFieldInterpolationCoodinates(points_coords, inv_points_coord_matrix); ++inv_quad_coord_it; ++inv_points_mat_it; ++quad_coords_it; ++points_coords_it; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::interpolateStress(const ElementType type, Vector & result) { AKANTU_DEBUG_IN(); #define INTERPOLATE_ELEMENTAL_FIELD(type) \ interpolateElementalField(stress(type), \ result) \ AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(INTERPOLATE_ELEMENTAL_FIELD); #undef INTERPOLATE_ELEMENTAL_FIELD AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::interpolateElementalField(const Vector & field, Vector & result) { AKANTU_DEBUG_IN(); Vector & elem_fil = element_filter(type); UInt nb_element = elem_fil.getSize(); UInt nb_quad_per_element = model->getFEM().getNbQuadraturePoints(type); UInt size_inverse_coords = getSizeElementalFieldInterpolationCoodinates(); types::RMatrix coefficients(nb_quad_per_element, field.getNbComponent()); const Vector & interp_inv_coord = interpolation_inverse_coordinates(type); const Vector & interp_points_coord = interpolation_points_matrices(type); UInt nb_interpolation_points_per_elem = interp_points_coord.getNbComponent() / size_inverse_coords; Vector::const_iterator field_it = field.begin_reinterpret(nb_quad_per_element, field.getNbComponent(), nb_element); Vector::const_iterator interpolation_points_coordinates_it = interp_points_coord.begin(nb_interpolation_points_per_elem, size_inverse_coords); Vector::iterator result_it = result.begin_reinterpret(nb_interpolation_points_per_elem, field.getNbComponent(), nb_element); Vector::const_iterator inv_quad_coord_it = interp_inv_coord.begin(size_inverse_coords, size_inverse_coords); /// loop over the elements of the current material and element type for (UInt el = 0; el < nb_element; ++el, ++field_it, ++result_it, ++inv_quad_coord_it, ++interpolation_points_coordinates_it) { /** * matrix containing the inversion of the quadrature points' * coordinates */ const types::RMatrix & inv_quad_coord_matrix = *inv_quad_coord_it; /** * multiply it by the field values over quadrature points to get * the interpolation coefficients */ coefficients.mul(inv_quad_coord_matrix, *field_it); /// matrix containing the points' coordinates const types::RMatrix & coord = *interpolation_points_coordinates_it; /// multiply the coordinates matrix by the coefficients matrix and store the result (*result_it).mul(coord, coefficients); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ const Vector & Material::getVector(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) const { std::stringstream sstr; std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; sstr << id << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getVector(fvect_id); } catch(debug::Exception & e) { AKANTU_EXCEPTION("The material " << name << "(" < & Material::getVector(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) { std::stringstream sstr; std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; sstr << id << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getVector(fvect_id); } catch(debug::Exception & e) { AKANTU_EXCEPTION("The material " << name << "(" <(ByElementTypeVector & vect, UInt nb_component, bool temporary, ElementKind element_kind); template void Material::initInternalVector(ByElementTypeVector & vect, UInt nb_component, bool temporary, ElementKind element_kind); template void Material::initInternalVector(ByElementTypeVector & vect, UInt nb_component, bool temporary, ElementKind element_kind); template void Material::initInternalVector(ByElementTypeVector & vect, UInt nb_component, bool temporary, ElementKind element_kind); template void Material::resizeInternalVector(ByElementTypeVector & vect, ElementKind element_kind) const; template void Material::resizeInternalVector(ByElementTypeVector & vect, ElementKind element_kind) const; template void Material::resizeInternalVector(ByElementTypeVector & vect, ElementKind element_kind) const; template void Material::resizeInternalVector(ByElementTypeVector & vect, ElementKind element_kind) const; __END_AKANTU__ diff --git a/src/model/solid_mechanics/material_parameters.hh b/src/model/solid_mechanics/material_parameters.hh index f43696f62..839952852 100644 --- a/src/model/solid_mechanics/material_parameters.hh +++ b/src/model/solid_mechanics/material_parameters.hh @@ -1,143 +1,149 @@ /** * @file material_parameters.hh * * @author Nicolas Richart * * @date Thu Aug 09 21:06:34 2012 * * @brief class to handle the material parameters * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_PARAMETERS_HH__ #define __AKANTU_MATERIAL_PARAMETERS_HH__ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ enum ParamAccessType { _pat_internal = 0x0001, _pat_writable = 0x0010, _pat_readable = 0x0100, _pat_modifiable = 0x0110, //_pat_readable | _pat_writable, _pat_parsable = 0x1000, _pat_parsmod = 0x1110 }; + +ParamAccessType opertator|(const ParamAccessType & a, const ParamAccessType & b) { + ParamAccessType tmp = ParamAccessType(a | b); + return tmp; +} + template class MaterialParamTyped; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ class MaterialParam { public: MaterialParam(); MaterialParam(std::string name, std::string description, ParamAccessType param_type); virtual ~MaterialParam() {}; /* ------------------------------------------------------------------------ */ bool isInternal() const; bool isWritable() const; bool isReadable() const; bool isParsable() const; /* ------------------------------------------------------------------------ */ template void set(T & value); template const T & get() const; virtual void parseParam(std::string value); /* ------------------------------------------------------------------------ */ virtual void printself(std::ostream & stream) const; protected: template const MaterialParamTyped & getMaterialParamTyped() const; template MaterialParamTyped & getMaterialParamTyped(); private: std::string name; std::string description; ParamAccessType param_type; }; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template class MaterialParamTyped : public MaterialParam { public: MaterialParamTyped(std::string name, std::string description, ParamAccessType param_type, T & param); /* ------------------------------------------------------------------------ */ void setTyped(T & value); const T & getTyped() const; void parseParam(std::string value); virtual void printself(std::ostream & stream) const; private: T & param; }; /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ class MaterialParameters { public: ~MaterialParameters(); template void registerParam(std::string name, T & variable, ParamAccessType type, const std::string description = ""); template void registerParam(std::string name, T & variable, T default_value, ParamAccessType type, const std::string description = ""); /* ------------------------------------------------------------------------ */ template void set(std::string name, T value); template const T & get(std::string name) const; void parseParam(std::string name, std::string value); /* ------------------------------------------------------------------------ */ void printself(std::ostream & stream, int indent) const; private: std::map params; }; #include "material_parameters_tmpl.hh" __END_AKANTU__ #endif /* __AKANTU_MATERIAL_PARAMETERS_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.cc b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.cc index 863e323e4..5f8b3f57a 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.cc +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.cc @@ -1,267 +1,266 @@ /** * @file material_cohesive_exponential.cc * * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * * @date Mon Jul 09 14:13:56 2012 * * @brief Exponential irreversible cohesive law of mixed mode loading * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "material_cohesive_exponential.hh" #include "solid_mechanics_model.hh" #include "sparse_matrix.hh" #include "dof_synchronizer.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialCohesiveExponential::MaterialCohesiveExponential(SolidMechanicsModel & model, const ID & id) : MaterialCohesive(model,id) { AKANTU_DEBUG_IN(); this->registerParam("beta" , beta , 0. , _pat_parsable, "Beta parameter" ); this->registerParam("G_cI" , G_cI , 0. , _pat_parsable, "Mode I fracture energy" ); this->registerParam("G_cII" , G_cII , 0. , _pat_parsable, "Mode II fracture energy"); this->registerParam("kappa" , kappa , 0. , _pat_readable, "Kappa parameter" ); this->registerParam("delta_c", delta_c, 0. , _pat_readable, "Critical displacement" ); // this->initInternalVector(delta_max, 1, _ek_cohesive); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template MaterialCohesiveExponential::~MaterialCohesiveExponential() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::initMaterial() { AKANTU_DEBUG_IN(); MaterialCohesive::initMaterial(); kappa = G_cII / G_cI; delta_c = 2 * G_cI / sigma_c; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::resizeCohesiveVectors() { MaterialCohesive::resizeCohesiveVectors(); //this->resizeInternalVector(delta_max, _ek_cohesive); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::computeTraction(const Vector & normal, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators Vector::iterator traction_it = tractions(el_type, ghost_type).begin(spatial_dimension); Vector::iterator opening_it = opening(el_type, ghost_type).begin(spatial_dimension); Vector::const_iterator normal_it = normal.begin(spatial_dimension); Vector::iteratortraction_end = tractions(el_type, ghost_type).end(spatial_dimension); Vector::iteratordelta_max_it = delta_max(el_type, ghost_type).begin(); /// compute scalars Real beta2 = beta*beta; /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++opening_it, ++normal_it, ++delta_max_it) { /// compute normal and tangential opening vectors Real normal_opening_norm = opening_it->dot(*normal_it); types::Vector normal_opening(spatial_dimension); normal_opening = (*normal_it); normal_opening *= normal_opening_norm; types::Vector tangential_opening(spatial_dimension); tangential_opening = *opening_it; tangential_opening -= normal_opening; Real tangential_opening_norm = tangential_opening.norm(); /** * compute effective opening displacement * @f$ \delta = \sqrt{ * \frac{\beta^2}{\kappa^2} \Delta_t^2 + \Delta_n^2 } @f$ */ Real delta = tangential_opening_norm; delta *= delta * beta2; delta += normal_opening_norm * normal_opening_norm; delta = sqrt(delta); /// full damage case if (std::abs(delta) < Math::getTolerance()) { /// set traction to zero traction_it->clear(); } else { /// element not fully damaged /** * Compute traction loading @f$ \mathbf{T} = * e \sigma_c \frac{\delta}{\delta_c} e^{-\delta/ \delta_c}@f$ */ /** * Compute traction unloading @f$ \mathbf{T} = * \frac{t_{max}}{\delta_{max}} \delta @f$ */ *traction_it = tangential_opening; *traction_it *= beta2; *traction_it += normal_opening; /// crack opening case if (delta > *delta_max_it) { Real k = exp(1)*sigma_c*exp(- delta / delta_c)/delta_c; *traction_it *= k; /// update maximum displacement *delta_max_it = delta; } else { /// unloading-reloading case Real k = exp(1)*sigma_c*exp(- *delta_max_it / delta_c)/delta_c; *traction_it *= k; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveExponential::computeTangentTraction(const ElementType & el_type, Vector & tangent_matrix, const Vector & normal, GhostType ghost_type) { AKANTU_DEBUG_IN(); Vector::iterator tangent_it = tangent_matrix.begin(spatial_dimension, spatial_dimension); Vector::iterator tangent_end = tangent_matrix.end(spatial_dimension, spatial_dimension); Vector::const_iterator normal_it = normal.begin(spatial_dimension); Vector::iterator opening_it = opening(el_type, ghost_type).begin(spatial_dimension); Vector::iterator traction_it = tractions(el_type, ghost_type).begin(spatial_dimension); Vector::iteratordelta_max_it = delta_max(el_type, ghost_type).begin(); Real beta2 = beta*beta; /** * compute tangent matrix @f$ \frac{\partial \mathbf{t}} * {\partial \delta} = \hat{\mathbf{t}} \otimes * \frac{\partial (t/\delta)}{\partial \delta} * \frac{\hat{\mathbf{t}}}{\delta}+ \frac{t}{\delta} [ \beta^2 \mathbf{I} + * (1-\beta^2) (\mathbf{n} \otimes \mathbf{n})] @f$ **/ /** * In which @f$ * \frac{\partial(t/ \delta)}{\partial \delta} = * \left\{\begin{array} {l l} * -e \frac{\sigma_c}{\delta_c^2 }e^{-\delta / \delta_c} & \quad if * \delta \geq \delta_{max} \\ * 0 & \quad if \delta < \delta_{max}, \delta_n > 0 * \end{array}\right. @f$ **/ for (; tangent_it != tangent_end; ++tangent_it, ++normal_it, ++opening_it, ++traction_it) { Real normal_opening_norm = opening_it->dot(*normal_it); types::Vector normal_opening(spatial_dimension); normal_opening = (*normal_it); normal_opening *= normal_opening_norm; types::Vector tangential_opening(spatial_dimension); tangential_opening = *opening_it; tangential_opening -= normal_opening; + Real tangential_opening_norm = tangential_opening.norm(); Real delta = tangential_opening_norm; delta *= delta * beta2; delta += normal_opening_norm * normal_opening_norm; delta = sqrt(delta); types::RVector t_hat(tangential_opening); t_hat *= beta2; t_hat += normal_opening; types::RMatrix nn(spatial_dimension, spatial_dimension); nn.outerProduct(*normal_it, *normal_it); types::RMatrix I(spatial_dimension, spatial_dimension); I.eye(beta2); nn *= (1-beta2); I += nn; if(std::abs(delta) < Math::getTolerance()){ *tangent_it += I; *tangent_it *= exp(1)* sigma_c/delta_c; } else { Real traction_norm = traction_it->norm(); I *= traction_norm / delta; types::RVector t_hat_tmp (t_hat); Real temp_var = 0; if ((delta > *delta_max_it) || (std::abs(delta - *delta_max_it) < 1e-12)) { temp_var = -exp(1- delta/delta_c) * sigma_c/(delta_c * delta_c); } temp_var /= delta; t_hat_tmp *= temp_var; types::RMatrix t_var_t(spatial_dimension, spatial_dimension); t_var_t.outerProduct(t_hat, t_hat_tmp); *tangent_it += I; *tangent_it += t_var_t; } } AKANTU_DEBUG_OUT(); - } -/* -------------------------------------------------------------------------- */ INSTANSIATE_MATERIAL(MaterialCohesiveExponential); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.hh b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.hh index fd995f34b..45269917a 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.hh +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.hh @@ -1,133 +1,131 @@ /** * @file material_cohesive_exponential.hh * * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * * @date Mon Jul 09 14:13:56 2012 * * @brief Exponential irreversible cohesive law of mixed mode loading * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "material_cohesive.hh" #include "aka_common.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_COHESIVE_EXPONENTIAL_HH__ #define __AKANTU_MATERIAL_COHESIVE_EXPONENTIAL_HH__ /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /** * Cohesive material Exponential damage * * parameters in the material files : * - sigma_c : critical stress sigma_c (default: 0) * - beta : weighting parameter for sliding and normal opening (default: 0) * - G_cI : fracture energy for mode I (default: 0) * - G_cII : fracture energy for mode II (default: 0) */ template class MaterialCohesiveExponential : public MaterialCohesive { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialCohesiveExponential(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialCohesiveExponential(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize the material computed parameter virtual void initMaterial(); /// resize vectors for new cohesive elements virtual void resizeCohesiveVectors(); protected: /// constitutive law void computeTraction(const Vector & normal, ElementType el_type, GhostType ghost_type = _not_ghost); /// compute the tangent stiffness matrix for an element type void computeTangentTraction(const ElementType & el_type, Vector & tangent_matrix, const Vector & normal, GhostType ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: - /// critical stress - Real sigma_c; /// beta parameter Real beta; /// mode I fracture energy Real G_cI; /// mode II fracture energy Real G_cII; /// kappa parameter Real kappa; // /// maximum displacement // ByElementTypeReal delta_max; /// critical displacement Real delta_c; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_cohesive_exponential_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_MATERIAL_COHESIVE_EXPONENTIAL_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_extrinsic.cc b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_extrinsic.cc index 8f6fef221..7808d92a1 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_extrinsic.cc +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_extrinsic.cc @@ -1,319 +1,319 @@ /** * @file material_cohesive_linear_extrinsic.cc * * @author Marco Vocialta * * @date Tue May 08 13:01:18 2012 * * @brief Linear irreversible cohesive law of mixed mode loading with * random stress definition for extrinsic type * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "material_cohesive_linear_extrinsic.hh" #include "solid_mechanics_model_cohesive.hh" #include "sparse_matrix.hh" #include "dof_synchronizer.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialCohesiveLinearExtrinsic::MaterialCohesiveLinearExtrinsic(SolidMechanicsModel & model, const ID & id) : MaterialCohesive(model,id), sigma_c_eff("sigma_c_eff",id), delta_c("delta_c",id), normal_stress(spatial_dimension), tangential_stress(spatial_dimension) { AKANTU_DEBUG_IN(); this->registerParam("beta" , beta , 0. , _pat_parsable, "Beta parameter" ); this->registerParam("G_cI" , G_cI , 0. , _pat_parsable, "Mode I fracture energy" ); this->registerParam("G_cII" , G_cII , 0. , _pat_parsable, "Mode II fracture energy"); this->registerParam("penalty", penalty, 0. , _pat_parsable, "Penalty coefficient" ); this->registerParam("kappa" , kappa , 1. , _pat_readable, "Kappa parameter" ); initInternalVector(sigma_c_eff, 1, false, _ek_cohesive); initInternalVector( delta_c, 1, false, _ek_cohesive); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template MaterialCohesiveLinearExtrinsic::~MaterialCohesiveLinearExtrinsic() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearExtrinsic::initMaterial() { AKANTU_DEBUG_IN(); MaterialCohesive::initMaterial(); if (G_cII != 0) kappa = G_cII / G_cI; /// compute scalars beta2_kappa2 = beta * beta/kappa/kappa; beta2_kappa = beta * beta/kappa; if (beta == 0) beta2_inv = 0; else beta2_inv = 1./beta/beta; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearExtrinsic::resizeCohesiveVectors() { MaterialCohesive::resizeCohesiveVectors(); - resizeInternalVector(sigma_c_eff, _ek_cohesive); - resizeInternalVector(delta_c, _ek_cohesive); + resizeInternalVector(sigma_c_eff, _ek_cohesive); + resizeInternalVector(delta_c, _ek_cohesive); FEM & fem_cohesive = model->getFEM("CohesiveFEM"); const Mesh & mesh = fem_cohesive.getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); for(; it != last_type; ++it) { const Vector & elem_filter = element_filter(*it, _not_ghost); UInt nb_element = elem_filter.getSize(); if (nb_element == 0) continue; UInt nb_quadrature_points = fem_cohesive.getNbQuadraturePoints(*it, _not_ghost); UInt nb_element_old = nb_element - sigma_insertion.getSize() / nb_quadrature_points; Vector & sigma_c_eff_vec = sigma_c_eff(*it, _not_ghost); Vector & delta_c_vec = delta_c(*it, _not_ghost); for (UInt el = nb_element_old; el < nb_element; ++el) { for (UInt q = 0; q < nb_quadrature_points; ++q) { Real new_sigma = sigma_insertion((el - nb_element_old) * nb_quadrature_points+q); Real new_delta = 2 * G_cI / new_sigma; sigma_c_eff_vec(el * nb_quadrature_points + q) = new_sigma; delta_c_vec(el * nb_quadrature_points + q) = new_delta; } } } } /* -------------------------------------------------------------------------- */ template inline Real MaterialCohesiveLinearExtrinsic::computeEffectiveNorm(const types::RMatrix & stress, const types::RVector & normal, const types::RVector & tangent) { AKANTU_DEBUG_IN(); normal_stress.mul(stress, normal); tangential_stress.mul(stress, tangent); Real normal_contrib = normal_stress.dot(normal); Real tangent_contrib = tangential_stress.dot(tangent); normal_contrib = std::max(0., normal_contrib); AKANTU_DEBUG_OUT(); return std::sqrt(normal_contrib*normal_contrib + tangent_contrib*tangent_contrib*beta2_inv); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearExtrinsic::computeStressNorms(const Vector & facet_stress, Vector & stress_check) { AKANTU_DEBUG_IN(); sigma_insertion.resize(0); Vector & facets_check = model->getFacetsCheck(); ElementType type_facet = model->getFacetType(); UInt nb_quad_facet = model->getFEM("FacetsFEM").getNbQuadraturePoints(type_facet); UInt nb_facet = facets_check.getSize(); const Vector & tangents = model->getTangents(); const Vector & normals = model->getFEM("FacetsFEM").getNormalsOnQuadPoints(type_facet); Vector::iterator stress_check_it = stress_check.begin(nb_quad_facet); Vector::const_iterator normal_it = normals.begin(spatial_dimension); Vector::const_iterator tangent_it = tangents.begin(spatial_dimension); Vector::const_iterator facet_stress_it = facet_stress.begin(spatial_dimension, spatial_dimension); for (UInt f = 0; f < nb_facet; ++f, ++stress_check_it) { for (UInt q = 0; q < nb_quad_facet; ++q, ++normal_it, ++tangent_it) { for (UInt e = 0; e < 2; ++e, ++facet_stress_it) { if (facets_check(f) == true) { Real effective_norm = computeEffectiveNorm(*facet_stress_it, *normal_it, *tangent_it); (*stress_check_it)(q) = std::max((*stress_check_it)(q), effective_norm); } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearExtrinsic::computeTraction(const Vector & normal, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators Vector::iterator traction_it = tractions(el_type, ghost_type).begin(spatial_dimension); Vector::iterator opening_it = opening(el_type, ghost_type).begin(spatial_dimension); Vector::const_iterator normal_it = normal.begin(spatial_dimension); Vector::iteratortraction_end = tractions(el_type, ghost_type).end(spatial_dimension); Vector::iteratorsigma_c_it = sigma_c_eff(el_type, ghost_type).begin(); Vector::iteratordelta_max_it = delta_max(el_type, ghost_type).begin(); Vector::iteratordelta_c_it = delta_c(el_type, ghost_type).begin(); Vector::iteratordamage_it = damage(el_type, ghost_type).begin(); Real epsilon = std::numeric_limits::epsilon(); Real * memory_space = new Real[2*spatial_dimension]; types::Vector normal_opening(memory_space, spatial_dimension); types::Vector tangential_opening(memory_space + spatial_dimension, spatial_dimension); /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++opening_it, ++normal_it, ++sigma_c_it, ++delta_max_it, ++delta_c_it, ++damage_it) { /// compute normal and tangential opening vectors Real normal_opening_norm = opening_it->dot(*normal_it); normal_opening = (*normal_it); normal_opening *= normal_opening_norm; tangential_opening = *opening_it; tangential_opening -= normal_opening; Real tangential_opening_norm = tangential_opening.norm(); /** * compute effective opening displacement * @f$ \delta = \sqrt{ * \frac{\beta^2}{\kappa^2} \Delta_t^2 + \Delta_n^2 } @f$ */ Real delta = tangential_opening_norm; delta *= delta * beta2_kappa2; /// don't consider penetration contribution if (normal_opening_norm > 0) delta += normal_opening_norm * normal_opening_norm; delta = sqrt(delta); /// full damage case or zero displacement case if (delta >= *delta_c_it || delta <= epsilon) { /// set traction to zero (*traction_it).clear(); if (normal_opening_norm < 0) { normal_opening *= penalty; *traction_it += normal_opening; } else { *damage_it = delta >= *delta_c_it; *delta_max_it = *damage_it * (*delta_c_it); } } /// element not fully damaged else { /** * Compute traction @f$ \mathbf{T} = \left( * \frac{\beta^2}{\kappa} \Delta_t \mathbf{t} + \Delta_n * \mathbf{n} \right) \frac{\sigma_c}{\delta} \left( 1- * \frac{\delta}{\delta_c} \right)@f$ */ *traction_it = tangential_opening; *traction_it *= beta2_kappa; /// update maximum displacement *delta_max_it = std::max(*delta_max_it, delta); *damage_it = *delta_max_it / *delta_c_it; Real k = *sigma_c_it / *delta_max_it * (1. - *damage_it); *traction_it *= k; /// use penalty coefficient in case of penetration if (normal_opening_norm < 0) k = penalty; normal_opening *= k; *traction_it += normal_opening; } } delete [] memory_space; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANSIATE_MATERIAL(MaterialCohesiveLinearExtrinsic); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc index 09a315bdc..e7156d684 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc +++ b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc @@ -1,637 +1,640 @@ /** * @file material_cohesive.cc * * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * * @date Wed Feb 22 16:31:20 2012 * * @brief Specialization of the material class for cohesive elements * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_cohesive.hh" #include "solid_mechanics_model_cohesive.hh" #include "sparse_matrix.hh" #include "dof_synchronizer.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ MaterialCohesive::MaterialCohesive(SolidMechanicsModel & model, const ID & id) : Material(model,id), reversible_energy("reversible_energy", id), total_energy("total_energy", id), tractions_old("tractions (old)",id), opening_old("opening (old)",id), tractions("tractions",id), opening("opening",id), delta_max("delta max",id), damage("damage", id) { AKANTU_DEBUG_IN(); this->model = dynamic_cast(&model); - initInternalVector(reversible_energy, 1, false, _ek_cohesive); - initInternalVector( total_energy, 1, false, _ek_cohesive); - initInternalVector( tractions_old, spatial_dimension, false, _ek_cohesive); - initInternalVector( tractions, spatial_dimension, false, _ek_cohesive); - initInternalVector( opening_old, spatial_dimension, false, _ek_cohesive); - initInternalVector( opening, spatial_dimension, false, _ek_cohesive); - initInternalVector( delta_max, 1, false, _ek_cohesive); - initInternalVector( damage, 1, false, _ek_cohesive); - this->registerParam("sigma_c", sigma_c, 0. , _pat_parsable, "Critical stress"); - this->registerParam("rand_factor", rand, 0. , _pat_parsable, "Randomness factor"); + this->registerParam("sigma_c", sigma_c, 0. , ParamAccessType(_pat_parsable | _pat_readable), "Critical stress"); + this->registerParam("rand_factor", rand, 0. , ParamAccessType(_pat_parsable | _pat_readable), "Randomness factor"); this->registerParam("distribution", distribution, std::string("uniform"), _pat_parsable, "Distribution type"); - this->registerParam("lambda", lambda, 0. , _pat_parsable, "Weibull modulus"); - this->registerParam("m", m_scale, 1. , _pat_parsable, "Scale parameter"); + this->registerParam("lambda", lambda, 0. , _pat_parsable, "Weibull modulus"); + this->registerParam("m", m_scale, 1. , _pat_parsable, "Scale parameter"); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MaterialCohesive::~MaterialCohesive() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::initMaterial() { AKANTU_DEBUG_IN(); Material::initMaterial(); + fem_cohesive = &(model->getFEMClass("CohesiveFEM")); Mesh & mesh = fem_cohesive->getMesh(); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, gt, _ek_cohesive); for(; it != last_type; ++it) { ElementType type = *it; element_filter.alloc(0, 1, type); } } + initInternalVector(reversible_energy, 1, false, _ek_cohesive); + initInternalVector( total_energy, 1, false, _ek_cohesive); + initInternalVector( tractions_old, spatial_dimension, false, _ek_cohesive); + initInternalVector( tractions, spatial_dimension, false, _ek_cohesive); + initInternalVector( opening_old, spatial_dimension, false, _ek_cohesive); + initInternalVector( opening, spatial_dimension, false, _ek_cohesive); + initInternalVector( delta_max, 1, false, _ek_cohesive); + initInternalVector( damage, 1, false, _ek_cohesive); + resizeCohesiveVectors(); - AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::resizeCohesiveVectors() { resizeInternalVector(reversible_energy, _ek_cohesive); - resizeInternalVector(total_energy, _ek_cohesive); - resizeInternalVector(tractions_old, _ek_cohesive); - resizeInternalVector(tractions, _ek_cohesive); - resizeInternalVector(opening_old, _ek_cohesive); - resizeInternalVector(opening, _ek_cohesive); - resizeInternalVector(delta_max, _ek_cohesive); - resizeInternalVector(damage, _ek_cohesive); + resizeInternalVector(total_energy , _ek_cohesive); + resizeInternalVector(tractions_old , _ek_cohesive); + resizeInternalVector(tractions , _ek_cohesive); + resizeInternalVector(opening_old , _ek_cohesive); + resizeInternalVector(opening , _ek_cohesive); + resizeInternalVector(delta_max , _ek_cohesive); + resizeInternalVector(damage , _ek_cohesive); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::generateRandomDistribution(Vector & sigma_lim) { AKANTU_DEBUG_IN(); std::srand(time(NULL)); UInt nb_facet = sigma_lim.getSize(); if (distribution == "uniform") { for (UInt i = 0; i < nb_facet; ++i) sigma_lim(i) = sigma_c * (1 + std::rand()/(Real)RAND_MAX * rand); } else if (distribution == "weibull") { Real exponent = 1./m_scale; for (UInt i = 0; i < nb_facet; ++i) sigma_lim(i) = sigma_c + lambda * std::pow(-1.* std::log(std::rand()/(Real)RAND_MAX), exponent); } else { AKANTU_DEBUG_ERROR("Unknown random distribution type for sigma_c"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::checkInsertion(const Vector & facet_stress, Vector & facet_insertion) { AKANTU_DEBUG_IN(); Vector & facets_check = model->getFacetsCheck(); ElementType type_facet = model->getFacetType(); UInt nb_quad_facet = model->getFEM("FacetsFEM").getNbQuadraturePoints(type_facet); UInt nb_facet = facets_check.getSize(); Vector stress_check(nb_facet, nb_quad_facet); stress_check.clear(); computeStressNorms(facet_stress, stress_check); bool * facet_check_it = facets_check.storage(); Vector::iterator stress_check_it = stress_check.begin(nb_quad_facet); Real * sigma_limit_it = model->getSigmaLimit().storage(); for (UInt f = 0; f < nb_facet; ++f, ++facet_check_it, ++stress_check_it, ++sigma_limit_it) { if (*facet_check_it == true) { for (UInt q = 0; q < nb_quad_facet; ++q) { if ((*stress_check_it)(q) > *sigma_limit_it) { facet_insertion.push_back(f); for (UInt qs = 0; qs < nb_quad_facet; ++qs) sigma_insertion.push_back((*stress_check_it)(qs)); *facet_check_it = false; break; } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the residual by assembling @f$\int_{e} t_e N_e dS @f$ * * @param[in] displacements nodes displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void MaterialCohesive::updateResidual(GhostType ghost_type) { AKANTU_DEBUG_IN(); /// compute traction computeTraction(ghost_type); /// update and assemble residual assembleResidual(ghost_type); /// compute energies computeEnergies(); /// update old values Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type, _ek_cohesive); for(; it != last_type; ++it) { tractions_old(*it, ghost_type).copy(tractions(*it, ghost_type)); opening_old(*it, ghost_type).copy(opening(*it, ghost_type)); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::assembleResidual(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Vector & residual = const_cast &>(model->getResidual()); Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type, _ek_cohesive); for(; it != last_type; ++it) { const Vector & shapes = fem_cohesive->getShapes(*it, ghost_type); Vector & elem_filter = element_filter(*it, ghost_type); Vector & traction = tractions(*it, ghost_type); UInt size_of_shapes = shapes.getNbComponent(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); UInt nb_quadrature_points = fem_cohesive->getNbQuadraturePoints(*it, ghost_type); UInt nb_element = elem_filter.getSize(); /// compute @f$t_i N_a@f$ Real * shapes_val = shapes.storage(); UInt * elem_filter_val = elem_filter.storage(); Vector * shapes_filtered = new Vector(nb_element*nb_quadrature_points, size_of_shapes, "filtered shapes"); Real * shapes_filtered_val = shapes_filtered->values; for (UInt el = 0; el < nb_element; ++el) { shapes_val = shapes.storage() + elem_filter_val[el] * size_of_shapes * nb_quadrature_points; memcpy(shapes_filtered_val, shapes_val, size_of_shapes * nb_quadrature_points * sizeof(Real)); shapes_filtered_val += size_of_shapes * nb_quadrature_points; } shapes_filtered_val = shapes_filtered->values; // multiply traction by shapes Vector * traction_cpy = new Vector(traction); traction_cpy->extendComponentsInterlaced(size_of_shapes, spatial_dimension); Real * traction_cpy_val = traction_cpy->storage(); for (UInt el = 0; el < nb_element; ++el) { for (UInt q = 0; q < nb_quadrature_points; ++q) { for (UInt n = 0; n < size_of_shapes; ++n,++shapes_filtered_val) { for (UInt i = 0; i < spatial_dimension; ++i) { *traction_cpy_val++ *= *shapes_filtered_val; } } } } delete shapes_filtered; /** * compute @f$\int t \cdot N\, dS@f$ by @f$ \sum_q \mathbf{N}^t * \mathbf{t}_q \overline w_q J_q@f$ */ Vector * int_t_N = new Vector(nb_element, spatial_dimension*size_of_shapes, "int_t_N"); fem_cohesive->integrate(*traction_cpy, *int_t_N, spatial_dimension*size_of_shapes, *it, ghost_type, &elem_filter); delete traction_cpy; int_t_N->extendComponentsInterlaced(2, int_t_N->getNbComponent() ); Real * int_t_N_val = int_t_N->storage(); for (UInt el = 0; el < nb_element; ++el) { for (UInt n = 0; n < size_of_shapes*spatial_dimension; ++n) int_t_N_val[n] *= -1.; int_t_N_val += nb_nodes_per_element*spatial_dimension; } /// assemble model->getFEMBoundary().assembleVector(*int_t_N, residual, model->getDOFSynchronizer().getLocalDOFEquationNumbers(), residual.getNbComponent(), *it, ghost_type, &elem_filter, 1); delete int_t_N; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::assembleStiffnessMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); SparseMatrix & K = const_cast(model->getStiffnessMatrix()); Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type, _ek_cohesive); for(; it != last_type; ++it) { UInt nb_quadrature_points = fem_cohesive->getNbQuadraturePoints(*it, ghost_type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); const Vector & shapes = fem_cohesive->getShapes(*it, ghost_type); Vector & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); UInt size_of_shapes = shapes.getNbComponent(); UInt * elem_filter_it = elem_filter.storage(); Vector * shapes_filtered = new Vector(nb_element*nb_quadrature_points, size_of_shapes, "filtered shapes"); Vector::iterator shapes_filtered_it = shapes_filtered->begin_reinterpret(size_of_shapes, nb_quadrature_points, nb_element); Vector::const_iterator shapes_it = shapes.begin_reinterpret(size_of_shapes, nb_quadrature_points, mesh.getNbElement(*it, ghost_type)); - for (UInt el = 0; el < nb_element; ++el) { + for (UInt el = 0; el < nb_element; ++el, ++shapes_filtered_it) { *shapes_filtered_it = shapes_it[elem_filter_it[el]]; } /** * compute A matrix @f$ \mathbf{A} = \left[\begin{array}{c c c c c c c c c c c c} * 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 \\ * 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 \\ * 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 \\ * 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 \\ * 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 \\ * 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 * \end{array} \right]@f$ **/ // UInt size_of_A = spatial_dimension*size_of_shapes*spatial_dimension*nb_nodes_per_element; // Real * A = new Real[size_of_A]; // memset(A, 0, size_of_A*sizeof(Real)); types::RMatrix A(spatial_dimension*size_of_shapes, spatial_dimension*nb_nodes_per_element); for ( UInt i = 0; i < spatial_dimension*size_of_shapes; ++i) { - A(i, i); + A(i, i) = 1; A(i, i + spatial_dimension*size_of_shapes) = -1; } /// compute traction computeTraction(ghost_type); /// get the tangent matrix @f$\frac{\partial{(t/\delta)}}{\partial{\delta}} @f$ Vector * tangent_stiffness_matrix = new Vector(nb_element*nb_quadrature_points, spatial_dimension* spatial_dimension, "tangent_stiffness_matrix"); - Vector * normal = new Vector(nb_element * nb_quadrature_points, - spatial_dimension, - "normal"); - computeNormal(model->getDisplacement(), *normal, *it, ghost_type); + // Vector * normal = new Vector(nb_element * nb_quadrature_points, spatial_dimension, "normal"); + Vector normal( nb_quadrature_points, spatial_dimension, "normal"); + + + computeNormal(model->getCurrentPosition(), normal, *it, ghost_type); tangent_stiffness_matrix->clear(); - computeTangentTraction(*it, *tangent_stiffness_matrix, *normal, ghost_type); + computeTangentTraction(*it, *tangent_stiffness_matrix, normal, ghost_type); - delete normal; + // delete normal; UInt size_at_nt_d_n_a = spatial_dimension*nb_nodes_per_element*spatial_dimension*nb_nodes_per_element; Vector * at_nt_d_n_a = new Vector (nb_element*nb_quadrature_points, size_at_nt_d_n_a, "A^t*N^t*D*N*A"); Vector::iterator > shapes_filt_it = shapes_filtered->begin(size_of_shapes); Vector::iterator D_it = tangent_stiffness_matrix->begin(spatial_dimension, spatial_dimension); Vector::iterator At_Nt_D_N_A_it = at_nt_d_n_a->begin(spatial_dimension * nb_nodes_per_element, spatial_dimension * nb_nodes_per_element); Vector::iterator At_Nt_D_N_A_end = at_nt_d_n_a->end (spatial_dimension * nb_nodes_per_element, spatial_dimension * nb_nodes_per_element); types::RMatrix N (spatial_dimension, spatial_dimension * size_of_shapes); types::RMatrix N_A (spatial_dimension, spatial_dimension * nb_nodes_per_element); types::RMatrix D_N_A(spatial_dimension, spatial_dimension * nb_nodes_per_element); for(; At_Nt_D_N_A_it != At_Nt_D_N_A_end; ++At_Nt_D_N_A_it, ++D_it, ++shapes_filt_it) { types::RMatrix & D = *D_it; types::RMatrix & At_Nt_D_N_A = *At_Nt_D_N_A_it; types::Vector & shapes_fil = *shapes_filt_it; N.clear(); /** * store the shapes in voigt notations matrix @f$\mathbf{N} = * \begin{array}{cccccc} N_0(\xi) & 0 & N_1(\xi) &0 & N_2(\xi) & 0 \\ * 0 & * N_0(\xi)& 0 &N_1(\xi)& 0 & N_2(\xi) \end{array} @f$ **/ for (UInt i = 0; i < spatial_dimension ; ++i) for (UInt n = 0; n < size_of_shapes; ++n) - N(i, i + n) = shapes_fil(n); + N(i, i + spatial_dimension * n) = shapes_fil(n); /** * compute stiffness matrix @f$ \mathbf{K} = \delta \mathbf{U}^T * \int_{\Gamma_c} {\mathbf{P}^t \frac{\partial{\mathbf{t}}} {\partial{\delta}} * \mathbf{P} d\Gamma \Delta \mathbf{U}} @f$ **/ N_A.mul(N, A); D_N_A.mul(D, N_A); At_Nt_D_N_A.mul(N_A, D_N_A); - } + } delete tangent_stiffness_matrix; delete shapes_filtered; Vector * K_e = new Vector(nb_element, size_at_nt_d_n_a, "K_e"); fem_cohesive->integrate(*at_nt_d_n_a, *K_e, size_at_nt_d_n_a, *it, ghost_type, &elem_filter); delete at_nt_d_n_a; model->getFEM().assembleMatrix(*K_e, K, spatial_dimension, *it, ghost_type, &elem_filter); delete K_e; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- * * Compute traction from displacements * * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void MaterialCohesive::computeTraction(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type, _ek_cohesive); for(; it != last_type; ++it) { Vector & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_quadrature_points = nb_element*fem_cohesive->getNbQuadraturePoints(*it, ghost_type); Vector normal(nb_quadrature_points, spatial_dimension, "normal"); /// compute normals @f$\mathbf{n}@f$ computeNormal(model->getCurrentPosition(), normal, *it, ghost_type); /// compute openings @f$\mathbf{\delta}@f$ computeOpening(model->getDisplacement(), opening(*it, ghost_type), *it, ghost_type); /// compute traction @f$\mathbf{t}@f$ computeTraction(normal, *it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeNormal(const Vector & position, Vector & normal, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); #define COMPUTE_NORMAL(type) \ fem_cohesive->getShapeFunctions(). \ computeNormalsOnControlPoints(position, \ normal, \ ghost_type, \ &(element_filter(type, ghost_type))) AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_NORMAL); #undef COMPUTE_NORMAL AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeOpening(const Vector & displacement, Vector & opening, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); #define COMPUTE_OPENING(type) \ fem_cohesive->getShapeFunctions(). \ interpolateOnControlPoints(displacement, \ opening, \ spatial_dimension, \ ghost_type, \ &(element_filter(type, ghost_type))) AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_OPENING); #undef COMPUTE_OPENING AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeEnergies() { AKANTU_DEBUG_IN(); Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); Real * memory_space = new Real[2*spatial_dimension]; types::RVector b(memory_space, spatial_dimension); types::RVector h(memory_space + spatial_dimension, spatial_dimension); for(; it != last_type; ++it) { Vector::iterator erev = reversible_energy(*it, _not_ghost).begin(); Vector::iterator etot = total_energy(*it, _not_ghost).begin(); Vector::iterator traction_it = tractions(*it, _not_ghost).begin(spatial_dimension); Vector::iterator traction_old_it = tractions_old(*it, _not_ghost).begin(spatial_dimension); Vector::iterator opening_it = opening(*it, _not_ghost).begin(spatial_dimension); Vector::iterator opening_old_it = opening_old(*it, _not_ghost).begin(spatial_dimension); Vector::iteratortraction_end = tractions(*it, _not_ghost).end(spatial_dimension); /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++traction_old_it, ++opening_it, ++opening_old_it, ++erev, ++etot) { /// trapezoidal integration b = *opening_it; b -= *opening_old_it; h = *traction_old_it; h += *traction_it; *etot += .5 * b.dot(h); *erev = .5 * traction_it->dot(*opening_it); } } delete [] memory_space; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getReversibleEnergy() { AKANTU_DEBUG_IN(); Real erev = 0.; /// integrate the dissipated energy for each type of elements Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); for(; it != last_type; ++it) { erev += fem_cohesive->integrate(reversible_energy(*it, _not_ghost), *it, _not_ghost, &element_filter(*it, _not_ghost)); } AKANTU_DEBUG_OUT(); return erev; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getDissipatedEnergy() { AKANTU_DEBUG_IN(); Real edis = 0.; /// integrate the dissipated energy for each type of elements Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); for(; it != last_type; ++it) { Vector dissipated_energy(total_energy(*it, _not_ghost)); dissipated_energy -= reversible_energy(*it, _not_ghost); edis += fem_cohesive->integrate(dissipated_energy, *it, _not_ghost, &element_filter(*it, _not_ghost)); } AKANTU_DEBUG_OUT(); return edis; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getEnergy(std::string type) { AKANTU_DEBUG_IN(); if (type == "reversible") return getReversibleEnergy(); else if (type == "dissipated") return getDissipatedEnergy(); AKANTU_DEBUG_OUT(); return 0.; } __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_elastic.cc b/src/model/solid_mechanics/materials/material_elastic.cc index a01245f6b..0b2a81752 100644 --- a/src/model/solid_mechanics/materials/material_elastic.cc +++ b/src/model/solid_mechanics/materials/material_elastic.cc @@ -1,113 +1,113 @@ /** * @file material_elastic.cc * * @author Nicolas Richart * * @date Tue Jul 27 18:15:37 2010 * * @brief Specialization of the material class for the elastic material * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "material_elastic.hh" #include "solid_mechanics_model.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialElastic::MaterialElastic(SolidMechanicsModel & model, const ID & id) : Material(model, id) { AKANTU_DEBUG_IN(); this->registerParam("E" ,E , 0. , ParamAccessType(_pat_parsable | _pat_modifiable), "Young's modulus" ); this->registerParam("nu" ,nu , 0.5 , ParamAccessType(_pat_parsable | _pat_modifiable), "Poisson's ratio" ); - this->registerParam("Plane_Stress",plane_stress, false, ParamAccessType(_pat_parsable | _pat_readable) , "Is plane stress" ); + this->registerParam("Plane_Stress",plane_stress, false, _pat_parsmod, "Is plane stress" ); /// @todo Plane_Stress should not be possible to be modified after initMaterial (but before) this->registerParam("lambda" ,lambda , _pat_readable, "First Lamé coefficient" ); this->registerParam("mu" ,mu , _pat_readable, "Second Lamé coefficient"); this->registerParam("kapa" ,kpa , _pat_readable, "Bulk coefficient" ); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialElastic::initMaterial() { AKANTU_DEBUG_IN(); Material::initMaterial(); if (spatial_dimension == 1) nu = 0.; updateInternalParameters(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialElastic::updateInternalParameters() { lambda = nu * E / ((1 + nu) * (1 - 2*nu)); mu = E / (2 * (1 + nu)); if(plane_stress) lambda = nu * E / ((1 + nu)*(1 - nu)); kpa = lambda + 2./3. * mu; } /* -------------------------------------------------------------------------- */ template void MaterialElastic::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); computeStressOnQuad(grad_u, sigma); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialElastic::computeTangentModuli(__attribute__((unused)) const ElementType & el_type, Vector & tangent_matrix, __attribute__((unused)) GhostType ghost_type) { AKANTU_DEBUG_IN(); MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix); computeTangentModuliOnQuad(tangent); MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template Real MaterialElastic::getPushWaveSpeed() const { return sqrt((lambda + 2*mu)/rho); } /* -------------------------------------------------------------------------- */ template Real MaterialElastic::getShearWaveSpeed() const { return sqrt(mu/rho); } /* -------------------------------------------------------------------------- */ INSTANSIATE_MATERIAL(MaterialElastic); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.cc b/src/model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.cc index 0687ebbe7..3ce23994c 100644 --- a/src/model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.cc +++ b/src/model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.cc @@ -1,176 +1,176 @@ /** * @file material_standard_linear_solid_deviatoric.cc * * @author Vladislav Yastrebov * @author Nicolas Richart * @author David Simon Kammer * * @date Wed Feb 08 16:53:31 2012 * * @brief Material Visco-elastic * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "material_standard_linear_solid_deviatoric.hh" #include "solid_mechanics_model.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialStandardLinearSolidDeviatoric::MaterialStandardLinearSolidDeviatoric(SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialElastic(model, id), stress_dev("stress_dev", id), history_integral("history_integral", id) { AKANTU_DEBUG_IN(); this->registerParam("Eta", eta, 1., ParamAccessType(_pat_parsable | _pat_modifiable), "Viscosity"); this->registerParam("Ev", Ev, 1., ParamAccessType(_pat_parsable | _pat_modifiable), "Stiffness of the viscous element"); - this->registerParam("Einf", E_inf, 1., ParamAccessType(_pat_readable | _pat_modifiable), "Stiffness of the elastic element"); + this->registerParam("Einf", E_inf, 1., ParamAccessType(_pat_readable), "Stiffness of the elastic element"); UInt stress_size = spatial_dimension * spatial_dimension; this->initInternalVector(this->stress_dev, stress_size); this->initInternalVector(this->history_integral, stress_size); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialStandardLinearSolidDeviatoric::initMaterial() { AKANTU_DEBUG_IN(); updateInternalParameters(); MaterialElastic::initMaterial(); this->resizeInternalVector(this->stress_dev); this->resizeInternalVector(this->history_integral); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialStandardLinearSolidDeviatoric::updateInternalParameters() { MaterialElastic::updateInternalParameters(); E_inf = this->E - this->Ev; } /* -------------------------------------------------------------------------- */ template void MaterialStandardLinearSolidDeviatoric::setToSteadyState(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Vector & stress_dev_vect = stress_dev(el_type, ghost_type); Vector & history_int_vect = history_integral(el_type, ghost_type); Vector::iterator stress_d = stress_dev_vect.begin(spatial_dimension, spatial_dimension); Vector::iterator history_int = history_int_vect.begin(spatial_dimension, spatial_dimension); /// Loop on all quadrature points MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); types::RMatrix & dev_s = *stress_d; types::RMatrix & h = *history_int; /// Compute the first invariant of strain Real Theta = grad_u.trace(); for (UInt i = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j) { dev_s(i, j) = 2 * this->mu * (.5 * (grad_u(i,j) + grad_u(j,i)) - 1./3. * Theta *(i == j)); h(i, j) = 0.; } /// Save the deviator of stress ++stress_d; ++history_int; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialStandardLinearSolidDeviatoric::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real tau = 0.; // if(std::abs(Ev) > std::numeric_limits::epsilon()) tau = eta / Ev; Vector & stress_dev_vect = stress_dev(el_type, ghost_type); Vector & history_int_vect = history_integral(el_type, ghost_type); Vector::iterator stress_d = stress_dev_vect.begin(spatial_dimension, spatial_dimension); Vector::iterator history_int = history_int_vect.begin(spatial_dimension, spatial_dimension); types::RMatrix s(spatial_dimension, spatial_dimension); Real dt = this->model->getTimeStep(); Real exp_dt_tau = exp( -dt/tau ); Real exp_dt_tau_2 = exp( -.5*dt/tau ); /// Loop on all quadrature points MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); types::RMatrix & dev_s = *stress_d; types::RMatrix & h = *history_int; s.clear(); sigma.clear(); /// Compute the first invariant of strain Real Theta = grad_u.trace(); Real gamma_inf = E_inf / this->E; Real gamma_v = Ev / this->E; types::RMatrix U_rond_prim(spatial_dimension, spatial_dimension); U_rond_prim.eye(gamma_inf * this->kpa * Theta); for (UInt i = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j) { s(i, j) = 2 * this->mu * (.5 * (grad_u(i,j) + grad_u(j,i)) - 1./3. * Theta *(i == j)); h(i, j) = exp_dt_tau * h(i, j) + exp_dt_tau_2 * (s(i, j) - dev_s(i, j)); dev_s(i, j) = s(i, j); sigma(i, j) = U_rond_prim(i,j) + gamma_inf * s(i, j) + gamma_v * h(i, j); } /// Save the deviator of stress ++stress_d; ++history_int; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANSIATE_MATERIAL(MaterialStandardLinearSolidDeviatoric); __END_AKANTU__ diff --git a/test/test_mesh_utils/test_mesh_io/test_mesh_io_msh.cc b/test/test_mesh_utils/test_mesh_io/test_mesh_io_msh.cc index d4b7fdec7..2ad4e6f01 100644 --- a/test/test_mesh_utils/test_mesh_io/test_mesh_io_msh.cc +++ b/test/test_mesh_utils/test_mesh_io/test_mesh_io_msh.cc @@ -1,52 +1,54 @@ /** * @file test_mesh_io_msh.cc * * @author Nicolas Richart * * @date Fri Sep 03 15:56:15 2010 * * @brief unit test for the MeshIOMSH class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { akantu::initialize(argc, argv); akantu::MeshIOMSH mesh_io; akantu::Mesh mesh(3); mesh_io.read("./cube.msh", mesh); std::cout << mesh << std::endl; + mesh_io.write("./cube.out", mesh); + akantu::finalize(); return EXIT_SUCCESS; } diff --git a/test/test_mesh_utils/test_mesh_partitionate/CMakeLists.txt b/test/test_mesh_utils/test_mesh_partitionate/CMakeLists.txt index 4c17c640c..2fb34851f 100644 --- a/test/test_mesh_utils/test_mesh_partitionate/CMakeLists.txt +++ b/test/test_mesh_utils/test_mesh_partitionate/CMakeLists.txt @@ -1,39 +1,45 @@ #=============================================================================== # @file CMakeLists.txt # # @author Guillaume Anciaux # # @date Fri Sep 03 15:56:15 2010 # # @brief configuration for mesh partitioner tests # # @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 . # # @section DESCRIPTION # #=============================================================================== add_mesh(test_mesh_partitionate_mesh triangle.geo 2 2) register_test(test_mesh_partitionate_scotch SOURCES test_mesh_partitionate_scotch.cc DEPENDENCIES test_mesh_partitionate_mesh DIRECTORIES_TO_CREATE paraview ) + +register_test(test_mesh_partitionate_scotch_advanced + SOURCE test_mesh_partitionate_scotch_advanced.cc + FILES_TO_COPY squares_L.msh squares_H.msh + DIRECTORIES_TO_CREATE paraview) + diff --git a/test/test_mesh_utils/test_mesh_partitionate/squares.geo b/test/test_mesh_utils/test_mesh_partitionate/squares.geo new file mode 100644 index 000000000..a25dfaafa --- /dev/null +++ b/test/test_mesh_utils/test_mesh_partitionate/squares.geo @@ -0,0 +1,43 @@ +// Element size +//h = 0.5; +h = 0.005; + +// Dimension of square +H = 2; +L = 1; + +// ------------------------------------------ +// Geometry +// ------------------------------------------ + +// Points +Point(1) = {-L/2., 0, 0, h}; +Point(2) = {L/2., 0, 0, h}; +Point(3) = {L/2., H, 0, h}; +Point(4) = {-L/2., H, 0, h}; + +Point(5) = {-L/2., 0, 0, h}; +Point(6) = {-L/2., -H, 0, h}; +Point(7) = {L/2., -H, 0, h}; +Point(8) = {L/2., 0, 0, h}; + +// Lines +Line(1) = {4, 1}; +Line(2) = {1, 2}; +Line(3) = {2, 3}; +Line(4) = {3, 4}; + +Line(5) = {8, 5}; +Line(6) = {5, 6}; +Line(7) = {6, 7}; +Line(8) = {7, 8}; + +// Geometric and Physical Surface +Line Loop(1) = {1, 2, 3, 4}; +Plane Surface(1) = {1}; + +Line Loop(2) = {5, 6, 7, 8}; +Plane Surface(2) = {2}; + +Transfinite Surface "*"; +Recombine Surface "*"; diff --git a/test/test_mesh_utils/test_mesh_partitionate/squares_H.msh b/test/test_mesh_utils/test_mesh_partitionate/squares_H.msh new file mode 100644 index 000000000..c0fdce822 --- /dev/null +++ b/test/test_mesh_utils/test_mesh_partitionate/squares_H.msh @@ -0,0 +1,87 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$Nodes +30 +1 -0.5 0 0 +2 0.5 0 0 +3 0.5 2 0 +4 -0.5 2 0 +5 -0.5 0 0 +6 -0.5 -2 0 +7 0.5 -2 0 +8 0.5 0 0 +9 -0.5 1.500000000000694 0 +10 -0.5 1.000000000002776 0 +11 -0.5 0.5000000000020817 0 +12 -1.040834085586084e-12 0 0 +13 0.5 0.4999999999987997 0 +14 0.5 0.9999999999980433 0 +15 0.5 1.499999999998675 0 +16 1.040834085586084e-12 2 0 +17 1.040834085586084e-12 0 0 +18 -0.5 -0.4999999999987997 0 +19 -0.5 -0.9999999999980433 0 +20 -0.5 -1.499999999998675 0 +21 -1.040834085586084e-12 -2 0 +22 0.5 -1.500000000000694 0 +23 0.5 -1.000000000002776 0 +24 0.5 -0.5000000000020817 0 +25 5.204170427930421e-13 1.499999999999684 0 +26 0 1.000000000000409 0 +27 -5.204309205808499e-13 0.5000000000004408 0 +28 5.20389287218676e-13 -0.5000000000004408 0 +29 2.036447582271663e-24 -1.000000000000409 0 +30 -5.203892872160472e-13 -1.499999999999685 0 +$EndNodes +$Elements +48 +1 15 2 0 1 1 +2 15 2 0 2 2 +3 15 2 0 3 3 +4 15 2 0 4 4 +5 15 2 0 5 5 +6 15 2 0 6 6 +7 15 2 0 7 7 +8 15 2 0 8 8 +9 1 2 0 1 4 9 +10 1 2 0 1 9 10 +11 1 2 0 1 10 11 +12 1 2 0 1 11 1 +13 1 2 0 2 1 12 +14 1 2 0 2 12 2 +15 1 2 0 3 2 13 +16 1 2 0 3 13 14 +17 1 2 0 3 14 15 +18 1 2 0 3 15 3 +19 1 2 0 4 3 16 +20 1 2 0 4 16 4 +21 1 2 0 5 8 17 +22 1 2 0 5 17 5 +23 1 2 0 6 5 18 +24 1 2 0 6 18 19 +25 1 2 0 6 19 20 +26 1 2 0 6 20 6 +27 1 2 0 7 6 21 +28 1 2 0 7 21 7 +29 1 2 0 8 7 22 +30 1 2 0 8 22 23 +31 1 2 0 8 23 24 +32 1 2 0 8 24 8 +33 3 2 0 1 4 9 25 16 +34 3 2 0 1 16 25 15 3 +35 3 2 0 1 9 10 26 25 +36 3 2 0 1 25 26 14 15 +37 3 2 0 1 10 11 27 26 +38 3 2 0 1 26 27 13 14 +39 3 2 0 1 11 1 12 27 +40 3 2 0 1 27 12 2 13 +41 3 2 0 2 8 17 28 24 +42 3 2 0 2 24 28 29 23 +43 3 2 0 2 23 29 30 22 +44 3 2 0 2 22 30 21 7 +45 3 2 0 2 17 5 18 28 +46 3 2 0 2 28 18 19 29 +47 3 2 0 2 29 19 20 30 +48 3 2 0 2 30 20 6 21 +$EndElements diff --git a/test/test_mesh_utils/test_mesh_partitionate/squares_L.msh b/test/test_mesh_utils/test_mesh_partitionate/squares_L.msh new file mode 100644 index 000000000..8919a07eb --- /dev/null +++ b/test/test_mesh_utils/test_mesh_partitionate/squares_L.msh @@ -0,0 +1,143 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$Nodes +54 +1 -2 0 0 +2 2 0 0 +3 2 1 0 +4 -2 1 0 +5 -2 0 0 +6 -2 -1 0 +7 2 -1 0 +8 2 0 0 +9 -2 0.5000000000013878 0 +10 -1.500000000001388 0 0 +11 -1.000000000002776 0 0 +12 -0.5000000000034694 0 0 +13 -4.163336342344337e-12 0 0 +14 0.4999999999965308 0 0 +15 0.9999999999972244 0 0 +16 1.499999999998612 0 0 +17 2 0.4999999999990217 0 +18 1.500000000001388 1 0 +19 1.000000000002776 1 0 +20 0.5000000000034694 1 0 +21 4.163336342344337e-12 1 0 +22 -0.4999999999965308 1 0 +23 -0.9999999999972244 1 0 +24 -1.499999999998612 1 0 +25 1.500000000001388 0 0 +26 1.000000000002776 0 0 +27 0.5000000000034694 0 0 +28 4.163336342344337e-12 0 0 +29 -0.4999999999965308 0 0 +30 -0.9999999999972244 0 0 +31 -1.499999999998612 0 0 +32 -2 -0.4999999999990217 0 +33 -1.500000000001388 -1 0 +34 -1.000000000002776 -1 0 +35 -0.5000000000034694 -1 0 +36 -4.163336342344337e-12 -1 0 +37 0.4999999999965308 -1 0 +38 0.9999999999972244 -1 0 +39 1.499999999998612 -1 0 +40 2 -0.5000000000013878 0 +41 -1.5 0.500000000001092 0 +42 -1 0.5000000000007963 0 +43 -0.4999999999999998 0.5000000000005005 0 +44 0 0.5000000000002047 0 +45 0.5000000000000001 0.499999999999909 0 +46 1 0.4999999999996132 0 +47 1.5 0.4999999999993174 0 +48 1.5 -0.500000000001092 0 +49 1 -0.5000000000007963 0 +50 0.5000000000000001 -0.5000000000005005 0 +51 8.145790329086651e-24 -0.5000000000002047 0 +52 -0.5000000000000001 -0.499999999999909 0 +53 -1 -0.4999999999996132 0 +54 -1.5 -0.4999999999993174 0 +$EndNodes +$Elements +80 +1 15 2 0 1 1 +2 15 2 0 2 2 +3 15 2 0 3 3 +4 15 2 0 4 4 +5 15 2 0 5 5 +6 15 2 0 6 6 +7 15 2 0 7 7 +8 15 2 0 8 8 +9 1 2 0 1 4 9 +10 1 2 0 1 9 1 +11 1 2 0 2 1 10 +12 1 2 0 2 10 11 +13 1 2 0 2 11 12 +14 1 2 0 2 12 13 +15 1 2 0 2 13 14 +16 1 2 0 2 14 15 +17 1 2 0 2 15 16 +18 1 2 0 2 16 2 +19 1 2 0 3 2 17 +20 1 2 0 3 17 3 +21 1 2 0 4 3 18 +22 1 2 0 4 18 19 +23 1 2 0 4 19 20 +24 1 2 0 4 20 21 +25 1 2 0 4 21 22 +26 1 2 0 4 22 23 +27 1 2 0 4 23 24 +28 1 2 0 4 24 4 +29 1 2 0 5 8 25 +30 1 2 0 5 25 26 +31 1 2 0 5 26 27 +32 1 2 0 5 27 28 +33 1 2 0 5 28 29 +34 1 2 0 5 29 30 +35 1 2 0 5 30 31 +36 1 2 0 5 31 5 +37 1 2 0 6 5 32 +38 1 2 0 6 32 6 +39 1 2 0 7 6 33 +40 1 2 0 7 33 34 +41 1 2 0 7 34 35 +42 1 2 0 7 35 36 +43 1 2 0 7 36 37 +44 1 2 0 7 37 38 +45 1 2 0 7 38 39 +46 1 2 0 7 39 7 +47 1 2 0 8 7 40 +48 1 2 0 8 40 8 +49 3 2 0 1 4 9 41 24 +50 3 2 0 1 24 41 42 23 +51 3 2 0 1 23 42 43 22 +52 3 2 0 1 22 43 44 21 +53 3 2 0 1 21 44 45 20 +54 3 2 0 1 20 45 46 19 +55 3 2 0 1 19 46 47 18 +56 3 2 0 1 18 47 17 3 +57 3 2 0 1 9 1 10 41 +58 3 2 0 1 41 10 11 42 +59 3 2 0 1 42 11 12 43 +60 3 2 0 1 43 12 13 44 +61 3 2 0 1 44 13 14 45 +62 3 2 0 1 45 14 15 46 +63 3 2 0 1 46 15 16 47 +64 3 2 0 1 47 16 2 17 +65 3 2 0 2 8 25 48 40 +66 3 2 0 2 40 48 39 7 +67 3 2 0 2 25 26 49 48 +68 3 2 0 2 48 49 38 39 +69 3 2 0 2 26 27 50 49 +70 3 2 0 2 49 50 37 38 +71 3 2 0 2 27 28 51 50 +72 3 2 0 2 50 51 36 37 +73 3 2 0 2 28 29 52 51 +74 3 2 0 2 51 52 35 36 +75 3 2 0 2 29 30 53 52 +76 3 2 0 2 52 53 34 35 +77 3 2 0 2 30 31 54 53 +78 3 2 0 2 53 54 33 34 +79 3 2 0 2 31 5 32 54 +80 3 2 0 2 54 32 6 33 +$EndElements diff --git a/test/test_mesh_utils/test_mesh_partitionate/test_mesh_partitionate_scotch_advanced.cc b/test/test_mesh_utils/test_mesh_partitionate/test_mesh_partitionate_scotch_advanced.cc new file mode 100644 index 000000000..2bdc1b061 --- /dev/null +++ b/test/test_mesh_utils/test_mesh_partitionate/test_mesh_partitionate_scotch_advanced.cc @@ -0,0 +1,391 @@ +/** + * @file test_mesh_partitionate_scotch_advanced.cc + * @author Nicolas Richart + * @author David Kammer + * @date Tue Oct 16 09:20:24 2012 + * + * @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 . + * + */ + +/* -------------------------------------------------------------------------- */ + +/* -------------------------------------------------------------------------- */ +#include "aka_common.hh" +#include "aka_vector.hh" +#include "mesh.hh" +#include "fem.hh" +#include "mesh_io_msh.hh" +#include "mesh_utils.hh" +#include "mesh_partition_scotch.hh" +/* -------------------------------------------------------------------------- */ +#ifdef AKANTU_USE_IOHELPER +# include "io_helper.hh" +#endif //AKANTU_USE_IOHELPER + +akantu::ElementType type = akantu::_quadrangle_4; +akantu::UInt nb_quadrature_points = 4; +akantu::UInt nb_nodes_per_element = 4; + +/* -------------------------------------------------------------------------- */ +void getInterfaceNodePairs(akantu::Mesh & mesh, + akantu::Vector & pairs); +void getPartition(const akantu::Mesh & mesh, + const akantu::MeshPartition::EdgeLoadFunctor & fctr, + const akantu::Vector & pairs, + akantu::MeshPartition * partition, + akantu::Vector & parts); +void dumpParaview(akantu::Mesh & mesh, std::string name, + double * part_1, double * part_2); + + +/* -------------------------------------------------------------------------- */ +class DoNotCutInterfaceFunctor : public akantu::MeshPartition::EdgeLoadFunctor { +public: + DoNotCutInterfaceFunctor(const akantu::Mesh & mesh) : mesh(mesh) { + + } + + virtual inline akantu::Int operator()(const akantu::Element & el1, + const akantu::Element & el2) const { + + const akantu::Vector & conn_1 = this->mesh.getConnectivity(el1.type); + const akantu::Vector & conn_2 = this->mesh.getConnectivity(el2.type); + std::set nodes; + + // count number of nodes + akantu::UInt nb_npe_1 = this->mesh.getNbNodesPerElement(el1.type); + for (akantu::UInt i=0; imesh.getNbNodesPerElement(el2.type); + for (akantu::UInt i=0; imesh.getNbElement(el1.type), + this->mesh.getNbElement(el2.type)); + int max_nb_npe = std::max(nb_npe_1, nb_npe_2); + + // get barycenter of elements to put different weights in vert. and horiz. dir. + akantu::Real bc_1[2]; + mesh.getBarycenter(el1.element,el1.type,bc_1); + akantu::Real bc_2[2]; + mesh.getBarycenter(el2.element,el2.type,bc_2); + akantu::Real bc_diff[2]; + bc_diff[0] = std::fabs(bc_1[0] - bc_2[0]); + bc_diff[1] = std::fabs(bc_1[1] - bc_2[1]); + + // put weights according to criterion + int weight = 1; + if (nodes.size() == nb_npe_1 + nb_npe_2) + weight = max_nb_element*max_nb_npe*4; + else if (bc_diff[0] < bc_diff[1]) + weight = 5; + + return weight; + } + +protected: + const akantu::Mesh & mesh; +}; + + +/* -------------------------------------------------------------------------- */ +/* Main */ +/* -------------------------------------------------------------------------- */ +int main(int argc, char *argv[]) +{ + akantu::initialize(argc, argv); + akantu::debug::setDebugLevel(akantu::dblDump); + akantu::debug::setDebugLevel(akantu::dblWarning); + + int dim = 2; + akantu::MeshIOMSH mesh_io; + + /* ---------- check if node pairs are considered with mesh_L --------- */ + akantu::Mesh mesh_l(dim); + mesh_io.read("squares_L.msh", mesh_l); + + // get interface node pairs + akantu::Vector pairs_l(0,2); + akantu::Vector pairs_empty(0,2); + getInterfaceNodePairs(mesh_l,pairs_l); + + akantu::MeshPartition * partition = new akantu::MeshPartitionScotch(mesh_l, mesh_l.getSpatialDimension()); + + // make normal partition -> it should cut along the interface + akantu::Vector parts(0,nb_quadrature_points); + getPartition(mesh_l, + akantu::MeshPartition::ConstEdgeLoadFunctor(), + pairs_empty, + partition, + parts); + double * part = parts.storage(); + + // make partition with node pairs -> it should cut perpendicular to the interface + akantu::Vector parts_adv(0,nb_quadrature_points); + getPartition(mesh_l, + akantu::MeshPartition::ConstEdgeLoadFunctor(), + pairs_l, + partition, + parts_adv); + double * part_adv = parts_adv.storage(); + + // output to visualize +#ifdef AKANTU_USE_IOHELPER + dumpParaview(mesh_l, "test-scotch-partition-mesh-L", + part, part_adv); +#endif //AKANTU_USE_IOHELPER + + // check + unsigned int nb_element = mesh_l.getNbElement(type); + akantu::Real bb_center[2]; + mesh_l.getBarycenter(0,type,bb_center); + + // define solution for part + unsigned int top_p_v = 0; + unsigned int bot_p_v = 1; + if (!(bb_center[1] > 0 && parts(0,0) == top_p_v) && + !(bb_center[1] < 0 && parts(0,0) == bot_p_v)) { + top_p_v = 1; + bot_p_v = 0; + } + std::cout << "top part = " << top_p_v << " | bot part = " << bot_p_v << std::endl; + + // define solution for part_adv + unsigned int left_p_v = 0; + unsigned int right_p_v = 1; + if (!(bb_center[0] > 0 && parts_adv(0,0) == right_p_v) && + !(bb_center[0] < 0 && parts_adv(0,0) == left_p_v)) { + left_p_v = 1; + right_p_v = 0; + } + std::cout << "left part = " << left_p_v << " | right part = " << right_p_v << std::endl; + + // check + for (akantu::UInt i=1; i 0 && parts(i,0) == top_p_v) && + !(bb_center[1] < 0 && parts(i,0) == bot_p_v)) { + std::cerr << " ** ERROR with mesh-L without node pairs" << std::endl; + return EXIT_FAILURE; + } + if (!(bb_center[0] > 0 && parts_adv(i,0) == right_p_v) && + !(bb_center[0] < 0 && parts_adv(i,0) == left_p_v)) { + std::cerr << " ** ERROR with mesh-L with node pairs" << std::endl; + return EXIT_FAILURE; + } + } + + delete partition; + + + /* ---------- check if node pairs and functor are considered with mesh_H --------- */ + akantu::Mesh mesh_h(dim,"mesh_h",1); + mesh_io.read("squares_H.msh", mesh_h); + akantu::Vector pairs_h(0,2); + getInterfaceNodePairs(mesh_h,pairs_h); + + partition = new akantu::MeshPartitionScotch(mesh_h, mesh_h.getSpatialDimension()); + + // make normal partition -> it should cut along the interface + getPartition(mesh_h, + akantu::MeshPartition::ConstEdgeLoadFunctor(), + pairs_h, + partition, + parts); + part = parts.storage(); + + // make partition with node pairs -> it should cut perpendicular to the interface + DoNotCutInterfaceFunctor fctr = DoNotCutInterfaceFunctor(mesh_h); + getPartition(mesh_h, + fctr, + pairs_h, + partition, + parts_adv); + part_adv = parts_adv.storage(); + + // output to visualize +#ifdef AKANTU_USE_IOHELPER + dumpParaview(mesh_h, "test-scotch-partition-mesh-H", + part, part_adv); +#endif //AKANTU_USE_IOHELPER + + // check + nb_element = mesh_h.getNbElement(type); + mesh_h.getBarycenter(0,type,bb_center); + + // define solution for part + top_p_v = 0; + bot_p_v = 1; + if (!(bb_center[1] > 0 && parts(0,0) == top_p_v) && + !(bb_center[1] < 0 && parts(0,0) == bot_p_v)) { + top_p_v = 1; + bot_p_v = 0; + } + std::cout << "top part = " << top_p_v << " | bot part = " << bot_p_v << std::endl; + + // define solution for part_adv + left_p_v = 0; + right_p_v = 1; + if (!(bb_center[0] > 0 && parts_adv(0,0) == right_p_v) && + !(bb_center[0] < 0 && parts_adv(0,0) == left_p_v)) { + left_p_v = 1; + right_p_v = 0; + } + std::cout << "left part = " << left_p_v << " | right part = " << right_p_v << std::endl; + + // check + for (akantu::UInt i=1; i 0 && parts(i,0) == top_p_v) && + !(bb_center[1] < 0 && parts(i,0) == bot_p_v)) { + std::cerr << " ** ERROR with mesh-H with node pairs and without functor" << std::endl; + return EXIT_FAILURE; + } + if (!(bb_center[0] > 0 && parts_adv(i,0) == right_p_v) && + !(bb_center[0] < 0 && parts_adv(i,0) == left_p_v)) { + std::cerr << " ** ERROR with mesh-H with node pairs and with functor" << std::endl; + return EXIT_FAILURE; + } + } + + delete partition; + + akantu::finalize(); + return EXIT_SUCCESS; +} + +/* -------------------------------------------------------------------------- */ +void dumpParaview(akantu::Mesh & mesh, std::string name, + double * part_1, double * part_2) { + + unsigned int nb_nodes = mesh.getNbNodes(); + unsigned int nb_element = mesh.getNbElement(type); + + iohelper::DumperParaview dumper; + dumper.SetMode(iohelper::TEXT); + dumper.SetPoints(mesh.getNodes().values, mesh.getSpatialDimension(), + nb_nodes, name); + dumper.SetConnectivity((int*) mesh.getConnectivity(type).values, + iohelper::QUAD1, nb_element, iohelper::C_MODE); + + dumper.AddElemDataField(part_1, 1, "partition_1"); + dumper.AddElemDataField(part_2, 1, "partition_2"); + dumper.SetPrefix("paraview"); + dumper.Init(); + dumper.Dump(); +} + +/* -------------------------------------------------------------------------- */ +void getPartition(const akantu::Mesh & mesh, + const akantu::MeshPartition::EdgeLoadFunctor & fctr, + const akantu::Vector & pairs, + akantu::MeshPartition * partition, + akantu::Vector & parts) { + + // partition->partitionate(2, akantu::MeshPartition::ConstEdgeLoadFunctor(), pairs_l); + partition->partitionate(2, fctr, pairs); + + // get partition value for each element + unsigned int nb_element = mesh.getNbElement(type); + + parts.resize(nb_element); + // double * part = new double[nb_element*nb_quadrature_points]; + akantu::UInt * part_val = partition->getPartition(type).values; + for (unsigned int i = 0; i < nb_element; ++i) + for (unsigned int q = 0; q < nb_quadrature_points; ++q) + //part[i*nb_quadrature_points + q] = part_val[i]; + parts(i,q) = part_val[i]; +} + + +/* -------------------------------------------------------------------------- */ +void getInterfaceNodePairs(akantu::Mesh & mesh, + akantu::Vector & pairs) { + + // put correct number of surfaces (gmsh starts with 1 but we need 0) + akantu::Vector & surf = const_cast &>(mesh.getUIntData(akantu::_segment_2, "tag_1")); + akantu::Vector::iterator<> it = surf.begin(); + akantu::Vector::iterator<> end = surf.end(); + for (;it != end; ++it) --(*it); + + // set surface id + mesh.setSurfaceIDsFromIntData("tag_1"); + akantu::CSR all_surface_nodes; + akantu::MeshUtils::buildNodesPerSurface(mesh, all_surface_nodes); + + // find top interface nodes + akantu::Vector top_nodes(0); + int top = 1; + for (akantu::CSR::iterator node = all_surface_nodes.begin(top); + node != all_surface_nodes.end(top); + ++node) { + top_nodes.push_back(*node); + } + + // find bottom interface nodes + akantu::Vector bot_nodes(0); + int bot = 4; + for (akantu::CSR::iterator node = all_surface_nodes.begin(bot); + node != all_surface_nodes.end(bot); + ++node) { + bot_nodes.push_back(*node); + } + + // verify that there is the same number of top and bottom nodes + int nb_pairs = 0; + if (bot_nodes.getSize() != top_nodes.getSize()) { + std::cerr << " ** ERROR: Interface does not have the same number of top and bottom nodes" << std::endl; + } + else { + nb_pairs = top_nodes.getSize(); + } + + // make pairs according to their x coordinate + const akantu::Vector & coords = mesh.getNodes(); + int dir = 0; + for (int i=0; i::max(); + double bot_min_v = std::numeric_limits::max(); + for (akantu::UInt j=0; j # @author Marco Vocialta # # @date Mon Jul 09 14:13:56 2012 # # @brief configuration for intrinsic implicit cohesive elements # # @section LICENSE # # Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== register_test(test_cohesive_intrinsic_impl test_cohesive_intrinsic_impl.cc) #=============================================================================== configure_file( ${CMAKE_CURRENT_SOURCE_DIR}/material.dat ${CMAKE_CURRENT_BINARY_DIR}/material.dat COPYONLY ) + +file(COPY material.dat DESTINATION .) file(COPY simple.msh DESTINATION .) file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/paraview) file(COPY simple.msh DESTINATION .) \ No newline at end of file diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic_impl/stiffness_matrix.verified b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic_impl/stiffness_matrix.verified new file mode 100644 index 000000000..81770e0ca --- /dev/null +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic_impl/stiffness_matrix.verified @@ -0,0 +1,308 @@ +%%MatrixMarket matrix coordinate real symmetric +36 36 306 +1 1 4.5 +1 2 -1.85037170770859e-16 +1 5 1 +1 6 8.32667268468867e-17 +1 31 -2.12792746386488e-16 +1 32 0.5 +1 3 -4 +1 4 -2.77555756156289e-16 +1 11 -1.55431223447522e-15 +1 12 -5.55111512312579e-17 +1 9 3.02535774210355e-15 +1 10 -2 +2 2 4.5 +2 5 -0.5 +2 6 0.5 +2 31 0.5 +2 32 -6.47630097698008e-17 +2 3 2 +2 4 -2 +2 11 0 +2 12 -6.66133814775094e-16 +2 9 -2 +2 10 1.99840144432528e-15 +5 5 4.5 +5 6 -1.5 +5 31 0.5 +5 32 0 +5 3 -4 +5 4 2 +5 11 -2 +5 12 2.22044604925031e-16 +5 9 7.7715611723761e-16 +5 10 -2.22044604925031e-16 +6 6 4.5 +6 31 -0.5 +6 32 1 +6 3 -3.19189119579733e-16 +6 4 -2 +6 11 2 +6 12 -4 +6 9 -2.35922392732846e-16 +6 10 3.33066907387547e-16 +31 31 4.5 +31 32 2.32317978986721e-16 +31 3 -1.11022302462516e-16 +31 4 -4.93038065763132e-32 +31 11 -2 +31 12 2 +31 9 1.19348975147204e-15 +31 10 -2 +32 32 5.10406262854645 +32 3 0 +32 4 -2.22044604925031e-16 +32 11 0 +32 12 -4 +32 9 -2 +32 10 9.29811783123569e-16 +3 3 12 +3 4 -2 +3 11 1.22124532708767e-15 +3 12 -2 +3 9 -4 +3 10 2 +4 4 12 +4 11 -2 +4 12 -8.88178419700125e-16 +4 9 2 +4 10 -8 +11 11 12 +11 12 -2 +11 9 -8 +11 10 2 +12 12 12 +12 9 2 +12 10 -4 +9 9 24 +9 10 -4 +10 10 24 +1 13 0.5 +1 14 -0.5 +1 15 -4.16333634234434e-16 +1 16 -1.11022302462516e-16 +1 7 -2 +1 8 2 +2 13 -6.47630097698008e-17 +2 14 1 +2 15 7.40148683083438e-17 +2 16 -1.0547118733939e-15 +2 7 4.44089209850063e-16 +2 8 -4 +31 13 1 +31 14 -1.23259516440783e-32 +31 15 -4 +31 16 4.64635957973443e-16 +31 7 2.77555756156289e-17 +31 8 0 +32 13 -0.5 +32 14 0.197968685726773 +32 15 2 +32 16 -1.39593737145355 +32 7 -5.55111512312578e-17 +32 8 1.38777878078145e-17 +13 13 4.5 +13 14 -1.5 +13 9 0 +13 10 0 +13 15 -4 +13 16 2 +13 7 -2 +13 8 0 +14 14 5.10406262854645 +14 9 -1.11022302462516e-16 +14 10 4.44089209850063e-16 +14 15 -3.53613655510927e-16 +14 16 -1.39593737145355 +14 7 2 +14 8 -4 +9 15 -4 +9 16 2 +9 7 -8 +9 8 2 +10 15 2 +10 16 -8 +10 7 2 +10 8 -4 +15 15 12 +15 16 -2 +15 7 1.77635683940025e-15 +15 8 -2 +16 16 14.4162505141858 +16 7 -2 +16 8 2.22044604925031e-15 +7 7 12 +7 8 -2 +8 8 12 +33 33 4.5 +33 34 1.10000786939369e-16 +33 17 1 +33 18 1.20274161001059e-16 +33 29 -3.05311331771918e-16 +33 30 0.5 +33 35 -4 +33 36 -9.64236319054764e-16 +33 23 -1.22124532708767e-15 +33 24 -1.29526019539602e-16 +33 21 7.07767178198537e-16 +33 22 -2 +34 34 5.10406262854646 +34 17 -0.5 +34 18 0.197968685726773 +34 29 0.5 +34 30 -5.55111512312577e-17 +34 35 2 +34 36 -1.39593737145355 +34 23 5.55111512312578e-17 +34 24 -5.55111512312578e-16 +34 21 -2 +34 22 -2.4980018054066e-16 +17 17 4.5 +17 18 -1.5 +17 29 0.5 +17 30 0 +17 35 -4 +17 36 2 +17 23 -2 +17 24 1.66533453693773e-16 +17 21 3.33066907387547e-16 +17 22 -1.66533453693773e-16 +18 18 5.10406262854645 +18 29 -0.5 +18 30 1 +18 35 2.00957989624968e-16 +18 36 -1.39593737145355 +18 23 2 +18 24 -4 +18 21 -2.22044604925031e-16 +18 22 5.55111512312578e-16 +29 29 4.5 +29 30 0 +29 35 -1.11022302462516e-16 +29 36 1.66533453693773e-16 +29 23 -2 +29 24 2 +29 21 2.88657986402541e-15 +29 22 -2 +30 30 4.5 +30 35 0 +30 36 -2.22044604925031e-16 +30 23 0 +30 24 -4 +30 21 -2 +30 22 1.60982338570648e-15 +35 35 12 +35 36 -2 +35 23 7.7715611723761e-16 +35 24 -2 +35 21 -4 +35 22 2 +36 36 14.4162505141858 +36 23 -2 +36 24 0 +36 21 2 +36 22 -8 +23 23 12 +23 24 -2 +23 21 -8 +23 22 2 +24 24 12 +24 21 2 +24 22 -4 +21 21 24 +21 22 -4 +22 22 24 +33 25 0.5 +33 26 -0.5 +33 27 2.35922392732846e-16 +33 28 -3.88578058618805e-16 +33 19 -2 +33 20 2 +34 25 3.70074341541719e-17 +34 26 1 +34 27 -2.96059473233375e-16 +34 28 -8.32667268468867e-17 +34 19 -2.46519032881566e-32 +34 20 -4 +29 25 1 +29 26 0 +29 27 -4 +29 28 0 +29 19 3.33066907387547e-16 +29 20 0 +30 25 -0.5 +30 26 0.5 +30 27 2 +30 28 -2 +30 19 0 +30 20 1.66533453693773e-16 +25 25 4.5 +25 26 -1.5 +25 21 6.66133814775094e-16 +25 22 -2.22044604925031e-16 +25 27 -4 +25 28 2 +25 19 -2 +25 20 -1.11022302462516e-16 +26 26 4.5 +26 21 -5.55111512312578e-16 +26 22 1.33226762955019e-15 +26 27 5.55111512312578e-16 +26 28 -2 +26 19 2 +26 20 -4 +21 27 -4.00000000000001 +21 28 2 +21 19 -8 +21 20 2 +22 27 2 +22 28 -8.00000000000001 +22 19 2 +22 20 -4 +27 27 12 +27 28 -2 +27 19 8.88178419700125e-16 +27 20 -2 +28 28 12 +28 19 -2 +28 20 1.33226762955019e-15 +19 19 12 +19 20 -2 +20 20 12 +31 17 -1.19130347991335e-31 +31 18 -2.32317978986721e-16 +31 33 5.95651739956675e-32 +31 34 1.23259516440783e-32 +31 35 -1.19130347991335e-31 +31 36 -4.64635957973443e-16 +32 17 -2.32317978986721e-16 +32 18 -0.604062628546454 +32 33 1.23259516440783e-32 +32 34 0.302031314273227 +32 35 -4.64635957973443e-16 +32 36 -0.604062628546454 +13 17 5.95651739956675e-32 +13 18 -1.23259516440783e-32 +13 33 -1.19130347991335e-31 +13 34 2.32317978986721e-16 +13 35 -1.19130347991335e-31 +13 36 4.64635957973443e-16 +14 17 -1.23259516440783e-32 +14 18 0.302031314273227 +14 33 2.32317978986721e-16 +14 34 -0.604062628546455 +14 35 4.64635957973443e-16 +14 36 -0.604062628546455 +15 17 -1.19130347991335e-31 +15 18 -4.64635957973443e-16 +15 33 -1.19130347991335e-31 +15 34 4.64635957973443e-16 +15 35 -4.7652139196534e-31 +15 36 0 +16 17 -4.64635957973443e-16 +16 18 -0.604062628546454 +16 33 4.64635957973443e-16 +16 34 -0.604062628546454 +16 35 0 +16 36 -2.41625051418582 diff --git a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic_impl/test_cohesive_intrinsic_impl.cc b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic_impl/test_cohesive_intrinsic_impl.cc index 0ed884856..14cbfbca7 100644 --- a/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic_impl/test_cohesive_intrinsic_impl.cc +++ b/test/test_model/test_solid_mechanics_model/test_cohesive/test_cohesive_intrinsic_impl/test_cohesive_intrinsic_impl.cc @@ -1,181 +1,183 @@ /** * @file test_cohesive_intrinsic_impl.cc * * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * * @date Mon Jul 09 14:13:56 2012 * * @brief Test for cohesive elements * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io_msh.hh" #include "mesh_utils.hh" #include "solid_mechanics_model_cohesive.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; int main(int argc, char *argv[]) { initialize(argc, argv); debug::setDebugLevel(dblError); const UInt spatial_dimension = 2; const ElementType type = _triangle_6; Mesh mesh(spatial_dimension); MeshIOMSH mesh_io; mesh_io.read("simple.msh", mesh); SolidMechanicsModelCohesive model(mesh); /// model initialization model.initFull("material.dat", _static); const Mesh & mesh_facets = model.getMeshFacets(); const ElementType type_facet = mesh.getFacetElementType(type); UInt nb_facet = mesh_facets.getNbElement(type_facet); Vector facet_insertion; Real * bary_facet = new Real[spatial_dimension]; for (UInt f = 0; f < nb_facet; ++f) { mesh_facets.getBarycenter(f, type_facet, bary_facet); - if (bary_facet[1] == 1) facet_insertion.push_back(f); + if (bary_facet[1] == 1){ + facet_insertion.push_back(f);} } delete[] bary_facet; model.insertCohesiveElements(facet_insertion); /// boundary conditions Vector & boundary = model.getBoundary(); UInt nb_nodes = mesh.getNbNodes(); Vector &position = const_cast &> (mesh.getNodes()); Vector & displacement = model.getDisplacement(); for (UInt n = 0; n < nb_nodes; ++n) { if (std::abs(position(n,1))< Math::getTolerance()){ boundary(n, 1) = true; displacement(n,1) = 0.0; } if ((std::abs(position(n,0))< Math::getTolerance())&& (position(n,1)< 1.1)){ boundary(n, 0) = true; displacement(n,0) = 0.0; } if ((std::abs(position(n,0)-1)< Math::getTolerance())&&(std::abs(position(n,1)-1)< Math::getTolerance())){ boundary(n, 0) = true; displacement(n,0) = 0.0; } if (std::abs(position(n,1)-2)< Math::getTolerance()){ boundary(n, 1) = true; } } model.setBaseName("intrinsic_impl"); model.addDumpField("displacement"); model.addDumpField("mass" ); model.addDumpField("velocity" ); model.addDumpField("acceleration"); model.addDumpField("force" ); model.addDumpField("residual" ); model.addDumpField("damage" ); model.addDumpField("stress" ); model.addDumpField("strain" ); model.dump(); const MaterialCohesive & mat_coh = dynamic_cast< const MaterialCohesive &> (model.getMaterial(1)); const Vector & opening = mat_coh.getOpening(_cohesive_2d_6); //const Vector & traction = mat_coh.getTraction(_cohesive_2d_6); - //model.getStiffnessMatrix().saveMatrix("K.mtx"); model.updateResidual(); const Vector & residual = model.getResidual(); UInt max_step = 1000; Real increment = 3./max_step; Real error_tol = 10e-6; std::ofstream fout; fout.open("output"); /// Main loop for ( UInt nstep = 0; nstep < max_step; ++nstep){ Real norm = 10; UInt count = 0; for (UInt n = 0; n < nb_nodes; ++n) { if (std::abs(position(n,1)-2)< Math::getTolerance()){ displacement(n,1) += increment; } } do{ std::cout << "Iter : " << ++count << " - residual norm : " << norm << std::endl; model.assembleStiffnessMatrix(); - model.getStiffnessMatrix().saveMatrix("K.mtx"); + if ((nstep == 0)&&(count == 1)) + model.getStiffnessMatrix().saveMatrix("stiffness_matrix.lastout"); + model.solveStatic(); model.updateResidual(); } while(!model.testConvergenceResidual(1e-5, norm) && (count < 100)) ; std::cout << "Step : " << nstep << " - residual norm : " << norm << std::endl; model.dump(); Real resid = 0; for (UInt n = 0; n < nb_nodes; ++n) { if (std::abs(position(n,1)-2)< Math::getTolerance()){ resid += residual.values[spatial_dimension* n + 1]; } } Real analytical = exp(1) * std::abs(opening.values[3]) * exp (-std::abs(opening.values[3])/0.5)/0.5; - // the residual force is comparing with the theoretical value of the cohesive law + //the residual force is comparing with the theoretical value of the cohesive law error_tol = std::abs((- resid - analytical)/analytical); fout << nstep << " " << -resid << " " << analytical << " " << error_tol << std::endl; if (error_tol > 1e-4) return EXIT_FAILURE; - } + } fout.close(); // finalize(); return EXIT_SUCCESS; }