diff --git a/src/common/aka_event_handler.hh b/src/common/aka_event_handler.hh index a1ec21ca1..8fad6f935 100644 --- a/src/common/aka_event_handler.hh +++ b/src/common/aka_event_handler.hh @@ -1,84 +1,85 @@ /** * @file aka_event_handler.hh * @author Nicolas Richart * @date Wed Aug 22 17:04:14 2012 * * @brief Base of Event Handler classes * * @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_AKA_EVENT_HANDLER_HH__ #define __AKANTU_AKA_EVENT_HANDLER_HH__ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ template class EventHandlerManager { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: + virtual ~EventHandlerManager() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void registerEventHandler(EventHandler & event_handler) { event_handlers.insert(&event_handler); } void unregisterEventHandler(EventHandler & event_handler) { event_handlers.erase(&event_handler); } template void sendEvent(const Event & event) { typename std::set::iterator it = event_handlers.begin(); typename std::set::iterator end = event_handlers.end(); for(;it != end; ++it) (*it)->sendEvent(event); } /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// list of the event handlers std::set event_handlers; }; __END_AKANTU__ #endif /* __AKANTU_AKA_EVENT_HANDLER_HH__ */ diff --git a/src/fem/mesh.hh b/src/fem/mesh.hh index ef675f4a8..99d2074b8 100644 --- a/src/fem/mesh.hh +++ b/src/fem/mesh.hh @@ -1,551 +1,553 @@ /** * @file mesh.hh * @author Nicolas Richart * @date Wed Jun 16 11:53:53 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)); } 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 { bool res = ((lhs.kind < rhs.kind) || ((lhs.kind == rhs.kind) && ((lhs.ghost_type < rhs.ghost_type) || ((lhs.ghost_type == rhs.ghost_type) && ((lhs.type < rhs.type) || ((lhs.type == rhs.type) && (lhs.element < rhs.element))))))); return res; } }; 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 NewNodesEvent : public MeshEvent { }; class RemovedNodesEvent : public MeshEvent { }; class NewElementsEvent : public MeshEvent { }; class RemovedElementsEvent : public MeshEvent { }; /* -------------------------------------------------------------------------- */ class MeshEventHandler { +public: + virtual ~MeshEventHandler() {}; /* ------------------------------------------------------------------------ */ /* Internal code */ /* ------------------------------------------------------------------------ */ public: inline void sendEvent(const NewNodesEvent & event) { onNodesAdded (event.getList()); } inline void sendEvent(const RemovedNodesEvent & event) { onNodesRemoved(event.getList()); } inline void sendEvent(const NewElementsEvent & event) { onElementsAdded (event.getList()); } inline void sendEvent(const RemovedElementsEvent & event) { onElementsRemoved(event.getList()); } /* ------------------------------------------------------------------------ */ /* Interface */ /* ------------------------------------------------------------------------ */ public: - virtual void onNodesAdded (const Vector & nodes_list) { } - virtual void onNodesRemoved(const Vector & nodes_list) { } + virtual void onNodesAdded (__attribute__((unused)) const Vector & nodes_list) { } + virtual void onNodesRemoved(__attribute__((unused)) const Vector & nodes_list) { } - virtual void onElementsAdded (const Vector & elements_list) { } - virtual void onElementsRemoved(const Vector & elements_list) { } + virtual void onElementsAdded (__attribute__((unused)) const Vector & elements_list) { } + virtual void onElementsRemoved(__attribute__((unused)) const Vector & elements_list) { } }; /* -------------------------------------------------------------------------- */ /* 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 size, const bool & flag_nb_node_per_elem_multiply = false, ElementKind element_kind = _ek_regular) 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); /// convert a linearized element to an element inline Element linearizedToElement (UInt linearized_element); /// 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); /* ------------------------------------------------------------------------ */ /* 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); /// @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/model/solid_mechanics/material.hh b/src/model/solid_mechanics/material.hh index 346939288..bb18ebf86 100644 --- a/src/model/solid_mechanics/material.hh +++ b/src/model/solid_mechanics/material.hh @@ -1,539 +1,545 @@ /** * @file material.hh * @author Nicolas Richart * @date Fri Jul 23 09:06:29 2010 * * @brief Mother class for all materials * * @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_memory.hh" #include "parser.hh" #include "data_accessor.hh" #include "material_parameters.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_HH__ #define __AKANTU_MATERIAL_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class Model; class SolidMechanicsModel; class CommunicationBuffer; } __BEGIN_AKANTU__ /** * Interface of all materials * Prerequisites for a new material * - inherit from this class * - implement the following methods: * \code * virtual Real getStableTimeStep(Real h, const Element & element = ElementNull); * * virtual void computeStress(ElementType el_type, * GhostType ghost_type = _not_ghost); * * virtual void computeTangentStiffness(const ElementType & el_type, * Vector & tangent_matrix, * GhostType ghost_type = _not_ghost); * \endcode * */ class Material : protected Memory, public DataAccessor, public Parsable, public MeshEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Material(SolidMechanicsModel & model, const ID & id = ""); virtual ~Material(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// register a material parameter template void registerParam(std::string name, T & variable, T default_value, ParamAccessType type, std::string description = ""); template void registerParam(std::string name, T & variable, ParamAccessType type, std::string description = ""); template - void registerInternal(ByElementTypeVector & vect) { AKANTU_DEBUG_TO_IMPLEMENT(); }; + void registerInternal(__attribute__((unused)) ByElementTypeVector & vect) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// read parameter from file virtual bool parseParam(const std::string & key, const std::string & value, const ID & id); /// function called to update the internal parameters when the modifiable /// parameters are modified virtual void updateInternalParameters() {} /// initialize the material computed parameter virtual void initMaterial(); /// compute the residual for this material virtual void updateResidual(GhostType ghost_type = _not_ghost); void assembleResidual(GhostType ghost_type); /// compute the residual for this material virtual void computeAllStresses(GhostType ghost_type = _not_ghost); - virtual void computeAllNonLocalStresses(GhostType ghost_type = _not_ghost) {}; + virtual void computeAllNonLocalStresses(__attribute__((unused)) GhostType ghost_type = _not_ghost) {}; /// set material to steady state void setToSteadyState(GhostType ghost_type = _not_ghost); /// compute the stiffness matrix virtual void assembleStiffnessMatrix(GhostType ghost_type); /// compute the stable time step for an element of size h virtual Real getStableTimeStep(__attribute__((unused)) Real h, __attribute__((unused)) const Element & element = ElementNull) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute the p-wave speed in the material virtual Real getPushWaveSpeed() const { AKANTU_DEBUG_TO_IMPLEMENT(); }; /// compute the s-wave speed in the material virtual Real getShearWaveSpeed() const { AKANTU_DEBUG_TO_IMPLEMENT(); }; /// add an element to the local mesh filter inline UInt addElement(const ElementType & type, UInt element, const GhostType & ghost_type); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; /** * interpolate stress on given positions for each element by means * of a geometrical interpolation on quadrature points */ virtual void interpolateStress(const ElementType type, Vector & result); /** * function to initialize the elemental field interpolation * function by inverting the quadrature points' coordinates */ virtual void initElementalFieldInterpolation(ByElementTypeReal & interpolation_points_coordinates); protected: /// constitutive law virtual void computeStress(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// set the material to steady state (to be implemented for materials that need it) virtual void setToSteadyState(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) {} /// compute the tangent stiffness matrix virtual void computeTangentModuli(__attribute__((unused)) const ElementType & el_type, __attribute__((unused)) Vector & tangent_matrix, __attribute__((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute the potential energy virtual void computePotentialEnergy(ElementType el_type, GhostType ghost_type = _not_ghost); template void assembleStiffnessMatrix(const ElementType & type, GhostType ghost_type); /// transfer the B matrix to a Voigt notation B matrix template inline void transferBMatrixToSymVoigtBMatrix(const types::Matrix & B, types::Matrix & Bvoigt, UInt nb_nodes_per_element) const; inline UInt getTangentStiffnessVoigtSize(UInt spatial_dimension) const; /// compute the potential energy by element void computePotentialEnergyByElement(); /// compute the coordinates of the quadrature points void computeQuadraturePointsCoordinates(ByElementTypeReal & quadrature_points_coordinates, const GhostType & ghost_type) const; /// interpolate an elemental field on given points for each element template void interpolateElementalField(const Vector & field, Vector & result); /// template function to initialize the elemental field interpolation template void initElementalFieldInterpolation(const Vector & quad_coordinates, const Vector & interpolation_points_coordinates); /// build the coordinate matrix for the interpolation on elemental field template inline void buildElementalFieldInterpolationCoodinates(const types::Matrix & coordinates, types::Matrix & coordMatrix); /// get the size of the coordiante matrix used in the interpolation template inline UInt getSizeElementalFieldInterpolationCoodinates(); /// extract the field values corresponding to the quadrature points used for the interpolation template inline void extractElementalFieldForInterplation(const Vector & field, Vector & filtered_field); /* ------------------------------------------------------------------------ */ /* Function for all materials */ /* ------------------------------------------------------------------------ */ protected: /// compute the potential energy for on element inline void computePotentialEnergyOnQuad(types::Matrix & grad_u, types::Matrix & sigma, Real & epot); protected: /// allocate an internal vector template void initInternalVector(ByElementTypeVector & vect, UInt nb_component, ElementKind element_kind = _ek_regular); public: /// resize an internal vector template void resizeInternalVector(ByElementTypeVector & vect, ElementKind element_kind = _ek_regular) const; /* ------------------------------------------------------------------------ */ template inline void gradUToF(const types::Matrix & grad_u, types::Matrix & F); inline void rightCauchy(const types::Matrix & F, types::Matrix & C); inline void leftCauchy (const types::Matrix & F, types::Matrix & B); /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ public: virtual inline UInt getNbDataToPack (const Element & element, SynchronizationTag tag) const; virtual inline UInt getNbDataToUnpack(const Element & element, SynchronizationTag tag) const; virtual UInt getNbDataToPack(__attribute__((unused)) SynchronizationTag tag) const { return 0; } virtual UInt getNbDataToUnpack(__attribute__((unused)) SynchronizationTag tag) const { return 0; } virtual inline void packData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag) const; virtual void packData(__attribute__((unused)) CommunicationBuffer & buffer, __attribute__((unused)) const UInt index, __attribute__((unused)) SynchronizationTag tag) const { } virtual inline void unpackData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag); virtual void unpackData(__attribute__((unused)) CommunicationBuffer & buffer, __attribute__((unused)) const UInt index, __attribute__((unused)) SynchronizationTag tag) { } /* ------------------------------------------------------------------------ */ virtual inline void onElementsAdded(const Vector & element_list); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Model, *model, const SolidMechanicsModel &) AKANTU_GET_MACRO(ID, id, const ID &); AKANTU_GET_MACRO(Rho, rho, Real); AKANTU_SET_MACRO(Rho, rho, Real); Real getPotentialEnergy(); virtual Real getEnergy(std::string type); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementFilter, element_filter, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Strain, strain, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(PotentialEnergy, potential_energy, Real); bool isNonLocal() const { return is_non_local; } const Vector & getVector(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost) const; Vector & getVector(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost); template inline T getParam(const ID & param) const; template inline void setParam(const ID & param, T value); protected: bool isInit() const { return is_init; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id of the material ID id; /// material name std::string name; /// The model to witch the material belong SolidMechanicsModel * model; /// density : rho Real rho; /// stresses arrays ordered by element types ByElementTypeReal stress; /// strains arrays ordered by element types ByElementTypeReal strain; /// list of element handled by the material ByElementTypeUInt element_filter; /// is the vector for potential energy initialized // bool potential_energy_vector; /// potential energy by element ByElementTypeReal potential_energy; /// tell if using in non local mode or not bool is_non_local; /// spatial dimension UInt spatial_dimension; /// elemental field interpolation coordinates ByElementTypeReal interpolation_inverse_coordinates; /// elemental field interpolation points ByElementTypeReal interpolation_points_matrices; /// list of the paramters MaterialParameters params; private: /// boolean to know if the material has been initialized bool is_init; std::map internal_vectors_real; std::map internal_vectors_uint; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #if defined (AKANTU_INCLUDE_INLINE_IMPL) # include "material_inline_impl.cc" #endif /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Material & _this) { _this.printself(stream); return stream; } __END_AKANTU__ /* -------------------------------------------------------------------------- */ /* Auto loop */ /* -------------------------------------------------------------------------- */ #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type) \ Vector::iterator strain_it = \ this->strain(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ Vector::iterator strain_end = \ this->strain(el_type, ghost_type).end(spatial_dimension, \ spatial_dimension); \ Vector::iterator stress_it = \ this->stress(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ \ for(;strain_it != strain_end; ++strain_it, ++stress_it) { \ types::Matrix & __attribute__((unused)) grad_u = *strain_it; \ types::Matrix & __attribute__((unused)) sigma = *stress_it #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END \ } \ #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_mat) \ Vector::iterator strain_it = \ this->strain(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ Vector::iterator strain_end = \ this->strain(el_type, ghost_type).end(spatial_dimension, \ spatial_dimension); \ \ UInt tangent_size = \ this->getTangentStiffnessVoigtSize(spatial_dimension); \ Vector::iterator tangent_it = \ tangent_mat.begin(tangent_size, \ tangent_size); \ \ for(;strain_it != strain_end; ++strain_it, ++tangent_it) { \ types::Matrix & __attribute__((unused)) grad_u = *strain_it; \ types::Matrix & tangent = *tangent_it #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END \ } /* -------------------------------------------------------------------------- */ /* Material list */ /* -------------------------------------------------------------------------- */ #define AKANTU_CORE_MATERIAL_LIST \ ((2, (elastic , MaterialElastic ))) \ ((2, (elastic_orthotropic, MaterialElasticOrthotropic ))) \ ((2, (neohookean , MaterialNeohookean ))) \ ((2, (sls_deviatoric , MaterialStandardLinearSolidDeviatoric))) \ ((2, (ve_stiffness_prop , MaterialStiffnessProportional ))) #define AKANTU_COHESIVE_MATERIAL_LIST \ ((2, (cohesive_bilinear , MaterialCohesiveBilinear ))) \ ((2, (cohesive_linear , MaterialCohesiveLinear ))) \ ((2, (cohesive_linear_extrinsic, MaterialCohesiveLinearExtrinsic ))) \ ((2, (cohesive_linear_exponential_extrinsic, MaterialCohesiveLinearExponentialExtrinsic ))) \ ((2, (cohesive_exponential , MaterialCohesiveExponential ))) #define AKANTU_DAMAGE_MATERIAL_LIST \ ((2, (damage_linear , MaterialDamageLinear ))) \ ((2, (marigo , MaterialMarigo ))) \ ((2, (mazars , MaterialMazars ))) \ ((2, (vreepeerlings , MaterialVreePeerlings ))) #ifdef AKANTU_DAMAGE_NON_LOCAL #define AKANTU_MATERIAL_WEIGHT_FUNCTION_TMPL_LIST \ ((stress_wf, StressBasedWeightFunction )) \ ((damage_wf, DamagedWeightFunction )) \ ((remove_wf, RemoveDamagedWeightFunction)) \ ((base_wf, BaseWeightFunction )) + +#define AKANTU_MATERIAL_VREEPEERLINGS_WEIGHT_FUNCTION_TMPL_LIST \ + AKANTU_MATERIAL_WEIGHT_FUNCTION_TMPL_LIST \ + ((removed_damrate_wf, RemoveDamagedWithDamageRateWeightFunction)) + #define AKANTU_DAMAGE_NON_LOCAL_MATERIAL_LIST \ ((3, (marigo_non_local , MaterialMarigoNonLocal, \ AKANTU_MATERIAL_WEIGHT_FUNCTION_TMPL_LIST))) \ - ((2, (mazars_non_local , MaterialMazarsNonLocal))) \ + ((2, (mazars_non_local , MaterialMazarsNonLocal ))) \ ((3, (vreepeerlings_non_local, MaterialVreePeerlingsNonLocal, \ - AKANTU_MATERIAL_WEIGHT_FUNCTION_TMPL_LIST))) + AKANTU_MATERIAL_VREEPEERLINGS_WEIGHT_FUNCTION_TMPL_LIST))) + #else # define AKANTU_DAMAGE_NON_LOCAL_MATERIAL_LIST #endif #define AKANTU_MATERIAL_LIST \ AKANTU_CORE_MATERIAL_LIST \ AKANTU_COHESIVE_MATERIAL_LIST \ AKANTU_DAMAGE_MATERIAL_LIST \ AKANTU_DAMAGE_NON_LOCAL_MATERIAL_LIST #define INSTANSIATE_MATERIAL(mat_name) \ template class mat_name<1>; \ template class mat_name<2>; \ template class mat_name<3> #if defined(__INTEL_COMPILER) #pragma warning ( push ) /* warning #654: overloaded virtual function "akantu::Material::computeStress" is only partially overridden in class "akantu::Material*" */ #pragma warning ( disable : 654 ) #endif //defined(__INTEL_COMPILER) /* -------------------------------------------------------------------------- */ // elastic materials #include "material_elastic.hh" #include "material_neohookean.hh" #include "material_elastic_orthotropic.hh" // visco-elastic materials #include "material_stiffness_proportional.hh" #include "material_standard_linear_solid_deviatoric.hh" // damage materials #include "material_damage.hh" #include "material_marigo.hh" #include "material_mazars.hh" #include "material_damage_linear.hh" #include "material_vreepeerlings.hh" #if defined(AKANTU_DAMAGE_NON_LOCAL) # include "material_marigo_non_local.hh" # include "material_mazars_non_local.hh" # include "material_vreepeerlings_non_local.hh" #endif // cohesive materials #include "material_cohesive.hh" #include "material_cohesive_linear.hh" #include "material_cohesive_bilinear.hh" #include "material_cohesive_linear_extrinsic.hh" #include "material_cohesive_exponential.hh" #include "material_cohesive_linear_exponential_extrinsic.hh" #if defined(__INTEL_COMPILER) #pragma warning ( pop ) #endif //defined(__INTEL_COMPILER) #endif /* __AKANTU_MATERIAL_HH__ */ diff --git a/src/model/solid_mechanics/material_inline_impl.cc b/src/model/solid_mechanics/material_inline_impl.cc index 064d05f34..d93822dad 100644 --- a/src/model/solid_mechanics/material_inline_impl.cc +++ b/src/model/solid_mechanics/material_inline_impl.cc @@ -1,268 +1,268 @@ /** * @file material_inline_impl.cc * @author Nicolas Richart * @date Tue Jul 27 11:57:43 2010 * * @brief Implementation of the inline functions of the class 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 . * */ __END_AKANTU__ #include "solid_mechanics_model.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ inline UInt Material::addElement(const ElementType & type, UInt element, const GhostType & ghost_type) { element_filter(type, ghost_type).push_back(element); return element_filter(type, ghost_type).getSize()-1; } /* -------------------------------------------------------------------------- */ inline UInt Material::getTangentStiffnessVoigtSize(UInt dim) const { return (dim * (dim - 1) / 2 + dim); } /* -------------------------------------------------------------------------- */ template inline void Material::gradUToF(const types::Matrix & grad_u, types::Matrix & F) { UInt size_F = F.size(); AKANTU_DEBUG_ASSERT(F.size() >= grad_u.size() && grad_u.size() == dim, "The dimension of the tensor F should be greater or equal to the dimension of the tensor grad_u."); for (UInt i = 0; i < dim; ++i) for (UInt j = 0; j < dim; ++j) F(i, j) = grad_u(i, j); for (UInt i = 0; i < size_F; ++i) F(i, i) += 1; } /* -------------------------------------------------------------------------- */ inline void Material::rightCauchy(const types::Matrix & F, types::Matrix & C) { C.mul(F, F); } /* -------------------------------------------------------------------------- */ inline void Material::leftCauchy(const types::Matrix & F, types::Matrix & B) { B.mul(F, F); } /* -------------------------------------------------------------------------- */ inline void Material::computePotentialEnergyOnQuad(types::Matrix & grad_u, types::Matrix & sigma, Real & epot) { epot = 0.; for (UInt i = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j) epot += sigma(i, j) * grad_u(i, j); epot *= .5; } /* -------------------------------------------------------------------------- */ template inline void Material::transferBMatrixToSymVoigtBMatrix(const types::Matrix & B, types::Matrix & Bvoigt, UInt nb_nodes_per_element) const { Bvoigt.clear(); for (UInt i = 0; i < dim; ++i) for (UInt n = 0; n < nb_nodes_per_element; ++n) Bvoigt(i, i + n*dim) = B(n, i); if(dim == 2) { ///in 2D, fill the @f$ [\frac{\partial N_i}{\partial x}, \frac{\partial N_i}{\partial y}]@f$ row for (UInt n = 0; n < nb_nodes_per_element; ++n) { Bvoigt(2, 1 + n*2) = B(n, 0); Bvoigt(2, 0 + n*2) = B(n, 1); } } if(dim == 3) { for (UInt n = 0; n < nb_nodes_per_element; ++n) { Real dndx = B(n, 0); Real dndy = B(n, 1); Real dndz = B(n, 2); ///in 3D, fill the @f$ [0, \frac{\partial N_i}{\partial y}, \frac{N_i}{\partial z}]@f$ row Bvoigt(3, 1 + n*3) = dndz; Bvoigt(3, 2 + n*3) = dndy; ///in 3D, fill the @f$ [\frac{\partial N_i}{\partial x}, 0, \frac{N_i}{\partial z}]@f$ row Bvoigt(4, 0 + n*3) = dndz; Bvoigt(4, 2 + n*3) = dndx; ///in 3D, fill the @f$ [\frac{\partial N_i}{\partial x}, \frac{N_i}{\partial y}, 0]@f$ row Bvoigt(5, 0 + n*3) = dndy; Bvoigt(5, 1 + n*3) = dndx; } } } /* -------------------------------------------------------------------------- */ template inline void Material::buildElementalFieldInterpolationCoodinates(__attribute__((unused)) const types::Matrix & coordinates, __attribute__((unused)) types::Matrix & coordMatrix) { AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ template<> inline void Material::buildElementalFieldInterpolationCoodinates<_triangle_3>(const types::Matrix & coordinates, types::Matrix & coordMatrix) { for (UInt i = 0; i < coordinates.rows(); ++i) coordMatrix(i, 0) = 1; } /* -------------------------------------------------------------------------- */ template<> inline void Material::buildElementalFieldInterpolationCoodinates<_triangle_6>(const types::Matrix & coordinates, types::Matrix & coordMatrix) { UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(_triangle_6); for (UInt i = 0; i < coordinates.rows(); ++i) { coordMatrix(i, 0) = 1; for (UInt j = 1; j < nb_quadrature_points; ++j) coordMatrix(i, j) = coordinates(i, j-1); } } /* -------------------------------------------------------------------------- */ template inline UInt Material::getSizeElementalFieldInterpolationCoodinates() { return model->getFEM().getNbQuadraturePoints(type); } /* -------------------------------------------------------------------------- */ template inline void Material::extractElementalFieldForInterplation(const Vector & field, Vector & filtered_field) { filtered_field.copy(field); } /* -------------------------------------------------------------------------- */ template void Material::registerParam(std::string name, T & variable, T default_value, ParamAccessType type, std::string description) { AKANTU_DEBUG_IN(); params.registerParam(name, variable, default_value, type, description); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::registerParam(std::string name, T & variable, ParamAccessType type, std::string description) { AKANTU_DEBUG_IN(); params.registerParam(name, variable, type, description); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline UInt Material::getNbDataToPack(const Element & element, SynchronizationTag tag) const { UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(element.type); if(tag == _gst_smm_stress) return spatial_dimension * spatial_dimension * sizeof(Real) * nb_quadrature_points; return 0; } /* -------------------------------------------------------------------------- */ inline UInt Material::getNbDataToUnpack(const Element & element, SynchronizationTag tag) const { UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(element.type); if(tag == _gst_smm_stress) return spatial_dimension * spatial_dimension * sizeof(Real) * nb_quadrature_points; return 0; } /* -------------------------------------------------------------------------- */ inline void Material::packData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag) const { if(tag == _gst_smm_stress) { UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(element.type); Vector::const_iterator stress_it = stress(element.type, _not_ghost).begin(spatial_dimension, spatial_dimension); stress_it += element.element * nb_quadrature_points; for (UInt q = 0; q < nb_quadrature_points; ++q, ++stress_it) buffer << *stress_it; } } /* -------------------------------------------------------------------------- */ inline void Material::unpackData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag) { if(tag == _gst_smm_stress) { UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(element.type); Vector::iterator stress_it = stress(element.type, _ghost).begin(spatial_dimension, spatial_dimension); stress_it += element.element * nb_quadrature_points; for (UInt q = 0; q < nb_quadrature_points; ++q, ++stress_it) buffer >> *stress_it; } } /* -------------------------------------------------------------------------- */ template inline T Material::getParam(const ID & param) const { try { return params.get(param); } catch (...) { AKANTU_EXCEPTION("No parameter " << param << " in the material " << id); } } /* -------------------------------------------------------------------------- */ template inline void Material::setParam(const ID & param, T value) { try { params.set(param, value); } catch(...) { AKANTU_EXCEPTION("No parameter " << param << " in the material " << id); } updateInternalParameters(); } /* -------------------------------------------------------------------------- */ -inline void Material::onElementsAdded(const Vector & element_list) { +inline void Material::onElementsAdded(__attribute__((unused)) const Vector & element_list) { for (std::map::iterator it = internal_vectors_real.begin(); it != internal_vectors_real.end(); ++it) { resizeInternalVector(*(it->second)); } for (std::map::iterator it = internal_vectors_uint.begin(); it != internal_vectors_uint.end(); ++it) { resizeInternalVector(*(it->second)); } } diff --git a/src/model/solid_mechanics/materials/material_damage/material_damage_non_local.hh b/src/model/solid_mechanics/materials/material_damage/material_damage_non_local.hh index ebf1e913a..94f7b23d8 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_damage_non_local.hh +++ b/src/model/solid_mechanics/materials/material_damage/material_damage_non_local.hh @@ -1,110 +1,136 @@ /** * @file material_damage_non_local.hh * @author Nicolas Richart * @date Thu Aug 16 18:28:00 2012 * * @brief interface for non local damage 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_DAMAGE_NON_LOCAL_HH__ #define __AKANTU_MATERIAL_DAMAGE_NON_LOCAL_HH__ __BEGIN_AKANTU__ template class MaterialDamageLocal, template class WeightFunction = BaseWeightFunction> class MaterialDamageNonLocal : public MaterialDamageLocal, public MaterialNonLocal { public: typedef MaterialNonLocal MaterialNonLocalParent; typedef MaterialDamageLocal MaterialDamageParent; MaterialDamageNonLocal(SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialDamageParent(model, id), MaterialNonLocalParent(model, id) { }; /* ------------------------------------------------------------------------ */ virtual void initMaterial() { MaterialDamageParent::initMaterial(); MaterialNonLocalParent::initMaterial(); } protected: /* -------------------------------------------------------------------------- */ virtual void computeNonLocalStress(ElementType type, GhostType ghost_type = _not_ghost) = 0; /* ------------------------------------------------------------------------ */ - void computeNonLocalStress(GhostType ghost_type) { + void computeNonLocalStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); Mesh::type_iterator it = this->model->getFEM().getMesh().firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = this->model->getFEM().getMesh().lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { computeNonLocalStress(*it, ghost_type); } this->updateDissipatedEnergy(ghost_type); AKANTU_DEBUG_OUT(); } public: /* ------------------------------------------------------------------------ */ virtual bool parseParam(const std::string & key, const std::string & value, const ID & id) { return MaterialNonLocalParent::parseParam(key, value, id) || MaterialDamageParent::parseParam(key, value, id); } public: /* ------------------------------------------------------------------------ */ virtual inline UInt getNbDataToPack(const Element & element, SynchronizationTag tag) const { return MaterialNonLocalParent::getNbDataToPack(element, tag) + MaterialDamageParent::getNbDataToPack(element, tag); } virtual inline UInt getNbDataToUnpack(const Element & element, SynchronizationTag tag) const { return MaterialNonLocalParent::getNbDataToUnpack(element, tag) + MaterialDamageParent::getNbDataToUnpack(element, tag); } + virtual inline UInt getNbDataToPack(SynchronizationTag tag) const { + return MaterialNonLocalParent::getNbDataToPack(tag) + + MaterialDamageParent::getNbDataToPack(tag); + } + + virtual inline UInt getNbDataToUnpack(SynchronizationTag tag) const { + return MaterialNonLocalParent::getNbDataToUnpack(tag) + + MaterialDamageParent::getNbDataToUnpack(tag); + } + virtual inline void packData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag) const { MaterialNonLocalParent::packData(buffer, element, tag); MaterialDamageParent::packData(buffer, element, tag); } virtual inline void unpackData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag) { MaterialNonLocalParent::unpackData(buffer, element, tag); MaterialDamageParent::unpackData(buffer, element, tag); } + + virtual void packData(CommunicationBuffer & buffer, + const UInt index, + SynchronizationTag tag) const { + MaterialNonLocalParent::packData(buffer, index, tag); + MaterialDamageParent::packData(buffer, index, tag); + } + + virtual void unpackData(CommunicationBuffer & buffer, + const UInt index, + SynchronizationTag tag) { + MaterialNonLocalParent::unpackData(buffer, index, tag); + MaterialDamageParent::unpackData(buffer, index, tag); + } + + }; __END_AKANTU__ #endif /* __AKANTU_MATERIAL_DAMAGE_NON_LOCAL_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_damage/material_marigo.cc b/src/model/solid_mechanics/materials/material_damage/material_marigo.cc index bfdd72154..7b19a763e 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_marigo.cc +++ b/src/model/solid_mechanics/materials/material_damage/material_marigo.cc @@ -1,125 +1,125 @@ /** * @file material_marigo.cc * @author Nicolas Richart * @author Guillaume Anciaux * @author Marion Chambart * @date Tue Jul 27 11:53:52 2010 * * @brief Specialization of the material class for the marigo 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_marigo.hh" #include "solid_mechanics_model.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialMarigo::MaterialMarigo(SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialDamage(model, id), Yd_rand("Yd_rand",id), damage_in_y(false), yc_limit(false) { AKANTU_DEBUG_IN(); this->registerParam("Yd", Yd, 50., ParamAccessType(_pat_parsable | _pat_modifiable)); this->registerParam("Sd", Sd, 5000., ParamAccessType(_pat_parsable | _pat_modifiable)); this->registerParam("Yd_randomness", Yd_randomness, 0., _pat_parsable, "Randomness in Yd"); this->registerParam("epsilon_c", epsilon_c, 0., _pat_parsable, "Critical strain"); this->registerParam("Yc limit", yc_limit, false, _pat_internal, "As the material a critical Y"); this->registerParam("damage_in_y", damage_in_y, false, _pat_parsable, "Use threshold (1-D)Y"); this->initInternalVector(this->Yd_rand, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialMarigo::initMaterial() { AKANTU_DEBUG_IN(); MaterialDamage::initMaterial(); updateInternalParameters(); this->resizeInternalVector(this->Yd_rand); const Mesh & mesh = this->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) { UInt nb_element = this->element_filter(*it).getSize(); UInt nb_quad = this->model->getFEM().getNbQuadraturePoints(*it); Vector & Yd_rand_vec = Yd_rand(*it); for(UInt e = 0; e < nb_element; ++e) { Real rand_part = (2 * drand48()-1) * Yd_randomness * Yd; for(UInt q = 0; q < nb_quad; ++q) Yd_rand_vec(nb_quad*e+q,0) = Yd + rand_part; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialMarigo::updateInternalParameters() { MaterialDamage::updateInternalParameters(); Yc = .5 * epsilon_c * this->E * epsilon_c; - if(epsilon_c != 0.) yc_limit = true; + if(std::abs(epsilon_c) > std::numeric_limits::epsilon()) yc_limit = true; else yc_limit = false; } /* -------------------------------------------------------------------------- */ template void MaterialMarigo::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real * dam = this->damage(el_type, ghost_type).storage(); Real * Yd_q = Yd_rand(el_type, ghost_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); Real Y = 0; computeStressOnQuad(grad_u, sigma, *dam, Y, *Yd_q); ++dam; ++Yd_q; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANSIATE_MATERIAL(MaterialMarigo); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_damage/material_mazars_non_local.cc b/src/model/solid_mechanics/materials/material_damage/material_mazars_non_local.cc index 4d6480396..58d9b7870 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_mazars_non_local.cc +++ b/src/model/solid_mechanics/materials/material_damage/material_mazars_non_local.cc @@ -1,176 +1,176 @@ /** * @file material_mazars_non_local.cc * @author Nicolas Richart * @author Marion Chambart * @date Tue Jul 27 11:53:52 2010 * * @brief Specialization of the material class for the non-local mazars 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_mazars_non_local.hh" #include "solid_mechanics_model.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialMazarsNonLocal::MaterialMazarsNonLocal(SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialMazars(model, id), MaterialNonLocalParent(model, id), Ehat("epsilon_equ", id) { AKANTU_DEBUG_IN(); this->damage_in_compute_stress = false; this->is_non_local = true; this->initInternalVector(this->Ehat, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialMazarsNonLocal::initMaterial() { AKANTU_DEBUG_IN(); MaterialMazars::initMaterial(); MaterialNonLocalParent::initMaterial(); this->resizeInternalVector(this->Ehat); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialMazarsNonLocal::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real * damage = this->damage(el_type, ghost_type).storage(); Real * epsilon_equ = this->Ehat(el_type, ghost_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); MaterialMazars::computeStressOnQuad(grad_u, sigma, *damage, *epsilon_equ); ++damage; ++epsilon_equ; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template -void MaterialMazarsNonLocal::computeNonLocalStress(GhostType ghost_type) { +void MaterialMazarsNonLocal::computeNonLocalStresses(GhostType ghost_type) { AKANTU_DEBUG_IN(); ByElementTypeReal nl_var("Non local variable", this->id); this->initInternalVector(nl_var, 1); this->resizeInternalVector(nl_var); if(this->damage_in_compute_stress) this->weightedAvergageOnNeighbours(this->damage, nl_var, 1); else this->weightedAvergageOnNeighbours(this->Ehat, nl_var, 1); Mesh::type_iterator it = this->model->getFEM().getMesh().firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = this->model->getFEM().getMesh().lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { this->computeNonLocalStress(nl_var(*it, ghost_type), *it, ghost_type); } this->updateDissipatedEnergy(ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialMazarsNonLocal::computeNonLocalStress(Vector & non_loc_var, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real * damage; Real * epsilon_equ; if(this->damage_in_compute_stress){ damage = non_loc_var.storage(); epsilon_equ = this->Ehat(el_type, ghost_type).storage(); } else { damage = this->damage(el_type, ghost_type).storage(); epsilon_equ = non_loc_var.storage(); } MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); this->computeDamageAndStressOnQuad(grad_u, sigma, *damage, *epsilon_equ); ++damage; ++epsilon_equ; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template bool MaterialMazarsNonLocal::parseParam(const std::string & key, const std::string & value, const ID & id) { std::stringstream sstr(value); if(key == "average_on_damage") { sstr >> this->damage_in_compute_stress; } else { return MaterialNonLocalParent::parseParam(key, value, id) || MaterialMazars::parseParam(key, value, id); } return true; } /* -------------------------------------------------------------------------- */ template void MaterialMazarsNonLocal::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "MaterialMazarsNonLocal [" << std::endl; if(this->damage_in_compute_stress) stream << space << " + Average on damage" << std::endl; else stream << space << " + Average on Ehat" << std::endl; MaterialMazars::printself(stream, indent + 1); MaterialNonLocalParent::printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ INSTANSIATE_MATERIAL(MaterialMazarsNonLocal); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_damage/material_mazars_non_local.hh b/src/model/solid_mechanics/materials/material_damage/material_mazars_non_local.hh index 1cbf63891..8e6e79fac 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_mazars_non_local.hh +++ b/src/model/solid_mechanics/materials/material_damage/material_mazars_non_local.hh @@ -1,103 +1,103 @@ /** * @file material_mazars_non_local.hh * @author * @date Wed Aug 31 17:08:23 2011 * * @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 "material.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_MAZARS_NON_LOCAL_HH__ #define __AKANTU_MATERIAL_MAZARS_NON_LOCAL_HH__ __BEGIN_AKANTU__ /** * Material Mazars Non local * * parameters in the material files : */ template class MaterialMazarsNonLocal : public MaterialMazars, public MaterialNonLocal { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef MaterialNonLocal MaterialNonLocalParent; MaterialMazarsNonLocal(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialMazarsNonLocal() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void initMaterial(); virtual bool parseParam(const std::string & key, const std::string & value, const ID & id); /// constitutive law for all element of a type void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); /// constitutive law - virtual void computeNonLocalStress(GhostType ghost_type = _not_ghost); - virtual void computeNonLocalStress(Vector & Ehatnl, - ElementType el_type, - GhostType ghost_type = _not_ghost); + void computeNonLocalStresses(GhostType ghost_type = _not_ghost); + void computeNonLocalStress(Vector & Ehatnl, + ElementType el_type, + GhostType ghost_type = _not_ghost); inline Real getStableTimeStep(Real h, const Element & element) { return MaterialMazars::getStableTimeStep(h, element); }; /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// the ehat per quadrature points to perform the averaging ByElementTypeReal Ehat; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "material_mazars_non_local_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_MATERIAL_MAZARS_NON_LOCAL_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_damage/material_vreepeerlings.cc b/src/model/solid_mechanics/materials/material_damage/material_vreepeerlings.cc index 6228fb9b0..4f4ac94fe 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_vreepeerlings.cc +++ b/src/model/solid_mechanics/materials/material_damage/material_vreepeerlings.cc @@ -1,208 +1,207 @@ /** * @file material_vreepeerlings.cc * @author Cyprien Wolff * @date Fri Feb 17 14:00:00 2012 * * @brief Specialization of the material class for the VreePeerlings 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_vreepeerlings.hh" #include "solid_mechanics_model.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialVreePeerlings::MaterialVreePeerlings(SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialDamage(model, id), - Kapa("Kapa",id) { + Kapa("Kapa",id), + strain_rate_vreepeerlings("strain-rate-vreepeerlings", id), + critical_strain("critical-strain", id) + { AKANTU_DEBUG_IN(); - Kapa0 = 0.0001; - Alpha = 0.99; - Beta = 300.; - Kct = 1.; - Kapa0_randomness = 0.; + this->registerParam("Kapa0i" , Kapa0i , 0.0001, _pat_parsable); + this->registerParam("Kapa0" , Kapa0 , 0.0001, _pat_parsable); + this->registerParam("Alpha" , Alpha , 0.99 , _pat_parsable); + this->registerParam("Beta" , Beta , 300. , _pat_parsable); + this->registerParam("Kct" , Kct , 1. , _pat_parsable); + this->registerParam("Kapa0_randomness", Kapa0_randomness, 0. , _pat_parsable); + + firststep = true; + countforinitialstep= 0; this->initInternalVector(this->Kapa, 1); + this->initInternalVector(this->critical_strain, 1); + this->initInternalVector(this->strain_rate_vreepeerlings, spatial_dimension * spatial_dimension); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialVreePeerlings::initMaterial() { AKANTU_DEBUG_IN(); MaterialDamage::initMaterial(); this->resizeInternalVector(this->Kapa); + this->resizeInternalVector(this->critical_strain); + this->resizeInternalVector(this->strain_rate_vreepeerlings); const Mesh & mesh = this->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) { Vector ::iterator kapa_it = Kapa(*it).begin(); Vector ::iterator kapa_end = Kapa(*it).end(); for(; kapa_it != kapa_end; ++kapa_it) { - Real rand_part = (2 * drand48()-1) * Kapa0_randomness * Kapa0; - *kapa_it = Kapa0 + rand_part; + Real rand_part = (2 * drand48()-1) * Kapa0_randomness * Kapa0i; + *kapa_it = Kapa0i + rand_part; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialVreePeerlings::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real * dam = this->damage(el_type, ghost_type).storage(); Real * Kapaq = Kapa(el_type, ghost_type).storage(); + Real * crit_strain = critical_strain(el_type, ghost_type).storage(); + Real dt = this->model->getTimeStep(); + + Vector & elem_filter = this->element_filter(el_type, ghost_type); + Vector & velocity = this->model->getVelocity(); + Vector & strain_rate_vrplgs = this->strain_rate_vreepeerlings(el_type, ghost_type); + + this->model->getFEM().gradientOnQuadraturePoints(velocity, strain_rate_vrplgs, + spatial_dimension, + el_type, ghost_type, &elem_filter); + Vector::iterator strain_rate_vrplgs_it = + strain_rate_vrplgs.begin(spatial_dimension, spatial_dimension); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); Real Equistrain; - computeStressOnQuad(grad_u, sigma, *dam, Equistrain, *Kapaq); + Real Equistrain_rate; + types::Matrix & strain_rate = *strain_rate_vrplgs_it; + + computeStressOnQuad(grad_u, sigma, *dam, Equistrain, Equistrain_rate, *Kapaq, dt, strain_rate, *crit_strain); ++dam; ++Kapaq; + ++strain_rate_vrplgs_it; + ++crit_strain; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; - AKANTU_DEBUG_OUT(); -} + if(!this->is_non_local) this->updateDissipatedEnergy(ghost_type); -/* -------------------------------------------------------------------------- */ -template -bool MaterialVreePeerlings::parseParam(const std::string & key, - const std::string & value, - const ID & id) { - std::stringstream sstr(value); - if(key == "Kapa0") { sstr >> Kapa0; } - else if(key == "Alpha") { sstr >> Alpha; } - else if(key == "Beta") { sstr >> Beta; } - else if(key == "Kct") { sstr >> Kct; } - else if(key == "Kapa0_randomness") { sstr >> Kapa0_randomness; } - else { return MaterialDamage::parseParam(key, value, id); } - return true; -} - - -/* -------------------------------------------------------------------------- */ -template -void MaterialVreePeerlings::printself(std::ostream & stream, int indent) const { - std::string space; - for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); - - stream << space << "Material<_vreepeerlings> [" << std::endl; - stream << space << " + Kapa0 : " << Kapa0 << std::endl; - stream << space << " + Alpha : " << Alpha << std::endl; - stream << space << " + Beta : " << Beta << std::endl; - stream << space << " + Kct : " << Kct << std::endl; - stream << space << " + Kapa0 randomness : " << Kapa0_randomness << std::endl; - MaterialDamage::printself(stream, indent + 1); - stream << space << "]" << std::endl; + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template UInt MaterialVreePeerlings::getNbDataToPack(const Element & element, SynchronizationTag tag) const { AKANTU_DEBUG_IN(); UInt size = 0; if(tag == _gst_smm_init_mat) { UInt nb_quad = this->model->getFEM().getNbQuadraturePoints(element.type); size += sizeof(Real) + nb_quad; } size += MaterialDamage::getNbDataToPack(element, tag); AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ template UInt MaterialVreePeerlings::getNbDataToUnpack(const Element & element, SynchronizationTag tag) const { AKANTU_DEBUG_IN(); UInt size = 0; if(tag == _gst_smm_init_mat) { UInt nb_quad = this->model->getFEM().getNbQuadraturePoints(element.type); size += sizeof(Real) + nb_quad; } size += MaterialDamage::getNbDataToPack(element, tag); AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ template void MaterialVreePeerlings::packData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag) const { AKANTU_DEBUG_IN(); if(tag == _gst_smm_init_mat){ UInt nb_quad = this->model->getFEM().getNbQuadraturePoints(element.type); const Vector & kapa = Kapa(element.type, _not_ghost); for(UInt q = 0; q < nb_quad; ++q) buffer << kapa(element.element * nb_quad + q); } MaterialDamage::packData(buffer, element, tag); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialVreePeerlings::unpackData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag) { AKANTU_DEBUG_IN(); if(tag == _gst_smm_init_mat) { UInt nb_quad = this->model->getFEM().getNbQuadraturePoints(element.type); Vector & kapa = Kapa(element.type, _not_ghost); for(UInt q = 0; q < nb_quad; ++q) buffer >> kapa(element.element * nb_quad + q); } MaterialDamage::packData(buffer, element, tag); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANSIATE_MATERIAL(MaterialVreePeerlings); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_damage/material_vreepeerlings.hh b/src/model/solid_mechanics/materials/material_damage/material_vreepeerlings.hh index d66034a80..b8015c671 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_vreepeerlings.hh +++ b/src/model/solid_mechanics/materials/material_damage/material_vreepeerlings.hh @@ -1,155 +1,170 @@ /** * @file material_vreepeerlings.hh * @author Cyprien Wolff * @date Fri Feb 17 14:00:00 2012 * * @brief Specialization of the material class for the VreePeerlings 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 "aka_common.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_VREEPEERLINGS_HH__ #define __AKANTU_MATERIAL_VREEPEERLINGS_HH__ __BEGIN_AKANTU__ /** * Material vreepeerlings * * parameters in the material files : - * - Kapa0 : (default: 50) Initial threshold (of the equivalent strain) + * - Kapa0i : (default: 0.0001) Initial threshold (of the equivalent strain) for the initial step + * - Kapa0 : (default: 0.0001) Initial threshold (of the equivalent strain) * - Alpha : (default: 0.99) Fitting parameter (must be close to 1 to do tend to 0 the stress in the damaged element) * - Beta : (default: 300) This parameter determines the rate at which the damage grows * - Kct : (default: 1) Ratio between compressive and tensile strength * - Kapa0_randomness : (default:0) Kapa random internal variable */ template class MaterialVreePeerlings : public MaterialDamage { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialVreePeerlings(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialVreePeerlings() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void initMaterial(); - bool parseParam(const std::string & key, const std::string & value, - const ID & id); - /// constitutive law for all element of a type void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); - /// function to print the containt of the class - virtual void printself(std::ostream & stream, int indent = 0) const; - protected: /// constitutive law for a given quadrature point inline void computeStressOnQuad(types::Matrix & F, types::Matrix & sigma, Real & dam, + Real & Equistrain_rate, Real & Equistrain, - Real & Kapaq); + Real & Kapaq, + Real dt, + types::Matrix & strain_rate_vrpgls, + Real & crit_strain); inline void computeDamageAndStressOnQuad(types::Matrix & sigma, Real & dam, + Real & Equistrain_rate, Real & Equistrain, - Real & Kapaq); + Real & Kapaq, + Real dt, + Real & crit_strain); /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ public: virtual UInt getNbDataToPack(const Element & element, SynchronizationTag tag) const; virtual UInt getNbDataToUnpack(const Element & element, SynchronizationTag tag) const; virtual void packData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag) const; virtual void unpackData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag); - virtual UInt getNbDataToPack(__attribute__((unused)) SynchronizationTag tag) const { return 0; } - virtual UInt getNbDataToUnpack(__attribute__((unused)) SynchronizationTag tag) const { return 0; } - virtual void packData(__attribute__((unused)) CommunicationBuffer & buffer, - __attribute__((unused)) const UInt index, - __attribute__((unused)) SynchronizationTag tag) const { } - virtual void unpackData(__attribute__((unused)) CommunicationBuffer & buffer, - __attribute__((unused)) const UInt index, - __attribute__((unused)) SynchronizationTag tag) { } + virtual UInt getNbDataToPack(SynchronizationTag tag) const { return MaterialDamage::getNbDataToPack(tag); } + virtual UInt getNbDataToUnpack(SynchronizationTag tag) const { return MaterialDamage::getNbDataToUnpack(tag); } + virtual void packData(CommunicationBuffer & buffer, + const UInt index, + SynchronizationTag tag) const { MaterialDamage::packData(buffer, index, tag); } + virtual void unpackData(CommunicationBuffer & buffer, + const UInt index, + SynchronizationTag tag) { MaterialDamage::unpackData(buffer, index, tag); } /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: + /// Initial threshold (of the equivalent strain) (used in the initial step) + Real Kapa0i; /// Initial threshold (of the equivalent strain) Real Kapa0; /// Fitting parameter (must be close to 1 to do tend to 0 the stress in the damaged element) Real Alpha; /// This parameter determines the rate at which the damage grows Real Beta; /// Ratio between compressive and tensile strength Real Kct; /// randomness on Kapa0 Real Kapa0_randomness; /// Kapa random internal variable ByElementTypeReal Kapa; + /// strain_rate_vreepeerlings + ByElementTypeReal strain_rate_vreepeerlings; + + /// strain_critical_vreepeerlings + ByElementTypeReal critical_strain; + + /// Booleen to check the first step + bool firststep; + + /// counter for initial step + Int countforinitialstep; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_vreepeerlings_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_MATERIAL_VREEPEERLINGS_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_damage/material_vreepeerlings_inline_impl.cc b/src/model/solid_mechanics/materials/material_damage/material_vreepeerlings_inline_impl.cc index a53acda24..e60cb2cee 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_vreepeerlings_inline_impl.cc +++ b/src/model/solid_mechanics/materials/material_damage/material_vreepeerlings_inline_impl.cc @@ -1,93 +1,250 @@ /** * @file material_vreepeerlings_inline_impl.cc * @author Cyprien Wolff * @date Fri Feb 17 14:00:00 2012 * * @brief Specialization of the material class for the VreePeerlings 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 . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template inline void MaterialVreePeerlings::computeStressOnQuad(types::Matrix & grad_u, types::Matrix & sigma, Real & dam, Real & Equistrain, - Real & Kapaq) { + Real & Equistrain_rate, + Real & Kapaq, + Real dt, + types::Matrix & strain_rate_vrplgs, + Real & crit_strain) { Real I1=0.; Real J2=0.; + Real I1rate=0.; + Real J2rate=0.; if(this->plane_stress) { I1 = (grad_u(0,0) + grad_u(1,1))*(1 - 2*this->nu)/(1 - this->nu); + I1rate = (strain_rate_vrplgs(0,0) + strain_rate_vrplgs(1,1))*(1 - 2*this->nu)/(1 - this->nu); + Real tmp = this->nu/(this->nu - 1); tmp *= tmp; + J2 = .5*(grad_u(0,0)*grad_u(0,0) + grad_u(1,1)*grad_u(1,1) + tmp*(grad_u(0,0) + grad_u(1,1))*(grad_u(0,0) + grad_u(1,1)) + .5*(grad_u(0,1) + grad_u(1,0))*(grad_u(0,1) + grad_u(1,0))) - I1*I1/6.; + // J2rate = .5*(strain_rate_vrplgs(0,0)*strain_rate_vrplgs(0,0) + + // strain_rate_vrplgs(1,1)*strain_rate_vrplgs(1,1) + + // tmp*(strain_rate_vrplgs(0,0) + strain_rate_vrplgs(1,1))*(strain_rate_vrplgs(0,0) + strain_rate_vrplgs(1,1)) + + // .5*(strain_rate_vrplgs(0,1) + strain_rate_vrplgs(1,0))*(strain_rate_vrplgs(0,1) + strain_rate_vrplgs(1,0))) - + // I1rate*I1rate/6.; + + J2rate = (grad_u(0,0)*strain_rate_vrplgs(0,0) + + grad_u(1,1)*strain_rate_vrplgs(1,1) + + tmp*(grad_u(0,0) + grad_u(1,1))*(strain_rate_vrplgs(0,0) + strain_rate_vrplgs(1,1)) + + (grad_u(0,1)*strain_rate_vrplgs(1,0)) + (grad_u(0,1)*strain_rate_vrplgs(1,0))) - + I1*I1rate/3.; + } else { I1 = grad_u.trace(); + I1rate = strain_rate_vrplgs.trace(); + J2 = .5*(grad_u(0,0)*grad_u(0,0) + grad_u(1,1)*grad_u(1,1) + grad_u(2,2)*grad_u(2,2) + .5*(grad_u(0,1) + grad_u(1,0))*(grad_u(0,1) + grad_u(1,0)) + .5*(grad_u(1,2) + grad_u(2,1))*(grad_u(1,2) + grad_u(2,1)) + .5*(grad_u(2,0) + grad_u(0,2))*(grad_u(2,0) + grad_u(0,2))) - I1 * I1/6.; + + // J2rate = .5*(strain_rate_vrplgs(0,0)*strain_rate_vrplgs(0,0) + + // strain_rate_vrplgs(1,1)*strain_rate_vrplgs(1,1) + + // strain_rate_vrplgs(2,2)*strain_rate_vrplgs(2,2) + + // .5*(strain_rate_vrplgs(0,1) + strain_rate_vrplgs(1,0))*(strain_rate_vrplgs(0,1) + strain_rate_vrplgs(1,0)) + + // .5*(strain_rate_vrplgs(1,2) + strain_rate_vrplgs(2,1))*(strain_rate_vrplgs(1,2) + strain_rate_vrplgs(2,1)) + + // .5*(strain_rate_vrplgs(2,0) + strain_rate_vrplgs(0,2))*(strain_rate_vrplgs(2,0) + strain_rate_vrplgs(0,2))) - + // I1rate * I1rate/6.; + + J2rate = (grad_u(0,0)*strain_rate_vrplgs(0,0) + + grad_u(1,1)*strain_rate_vrplgs(1,1) + + grad_u(2,2)*strain_rate_vrplgs(2,2) + + (grad_u(0,1)*strain_rate_vrplgs(1,0)) + (grad_u(0,1)*strain_rate_vrplgs(1,0)) + + (grad_u(1,2)*strain_rate_vrplgs(2,1)) + (grad_u(1,2)*strain_rate_vrplgs(2,1)) + + (grad_u(2,0)*strain_rate_vrplgs(0,2)) + (grad_u(2,0)*strain_rate_vrplgs(0,2))) - + I1 * I1rate/6.; + } + + //Real tmp = (Kct - 1)*I1/(1 - 2*this->nu); + Real tmp = (Kct - 1)/(1 - 2*this->nu); + Real tmp2 = (12 * Kct)/((1 + this->nu)*(1 + this->nu)); + // Real tmp3 = (12)/((1 + this->nu)*(1 + this->nu)); + ///////////////////////////////////////////////////debut + //Real tmprate = (Kct - 1)*I1rate/(1 - 2*this->nu); + ///////////////////////////////////////////////////fin + // Real nu1 = (1 + this->nu); + /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////debut + + Real Equistrainratenew=tmp*I1rate / (2*Kct) + 1./(4*Kct) * (2*tmp*tmp*I1*I1rate + tmp2*J2rate) / sqrt(tmp*tmp*I1*I1 + tmp2*J2); + Equistrain = tmp * I1 / (2*Kct) + 1./(2*Kct) * sqrt(tmp*tmp*I1*I1 + tmp2*J2); + + if(Equistrainratenew*Equistrain>0){ + if(I1rate * std::abs(I1rate) > 0){ + Equistrain_rate = Equistrainratenew; + }else{ + Equistrain_rate =tmp * std::abs(I1rate) / (2*Kct) + 1./(4*Kct)*(2*tmp*tmp*I1*I1rate + tmp2*J2rate) / sqrt(tmp*tmp*I1*I1 + tmp2*J2); + } + //test + //Equistrain_rate = 1./(4.) * (tmp3*J2rate) / sqrt(tmp3*J2); } - Real tmp = (Kct - 1)*I1/(1 - 2*this->nu); - Real nu1 = (1 + this->nu); + // Real tmprate = (Kct - 1)*abs(I1rate)/(1 - 2*this->nu); + // Real Equistrainratenew = tmprate / (2*Kct) + 1./(2*Kct) * sqrt(tmprate * tmprate + (12 * Kct * J2rate)/(nu1 * nu1)); + // if((I1*abs(I1)>0. && I1rate*abs(I1rate)>0.) || (I1*abs(I1)<0. && I1rate*abs(I1rate)<0.)){ + // Equistrain_rate=Equistrainratenew; + // } + /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////fin - Equistrain = tmp / (2*Kct) + 1./(2*Kct) * sqrt(tmp * tmp + (12 * Kct * J2)/(nu1 * nu1)); + // Equistrain = tmp / (2*Kct) + 1./(2*Kct) * sqrt(tmp * tmp + (12 * Kct * J2)/(nu1 * nu1)); + // Equistrain = tmp * I1 / (2*Kct) + 1./(2*Kct) * sqrt(tmp*tmp * I1*I1 + tmp2 * J2); + + /////////////////////////////////////////////////////////////////////////////////////////////////////////////debut + // Equistrain_rate = tmprate / (2*Kct) + 1./(2*Kct) * sqrt(tmprate * tmprate + (12 * Kct * J2rate)/(nu1 * nu1)); + /////////////////////////////////////////////////////////////////////////////////////////////////////////////fin MaterialElastic::computeStressOnQuad(grad_u, sigma); if(!this->is_non_local) { - computeDamageAndStressOnQuad(sigma, dam, Equistrain, Kapaq); + computeDamageAndStressOnQuad(sigma, dam, Equistrain, Equistrain_rate, Kapaq, dt, crit_strain); } } /* -------------------------------------------------------------------------- */ template inline void MaterialVreePeerlings::computeDamageAndStressOnQuad(types::Matrix & sigma, Real & dam, Real & Equistrain, - Real & Kapaq) { - Real Fd = Equistrain - Kapaq; - if (Fd > 0) { - Kapaq = std::max(Equistrain,Kapaq); - //dam = 1. - Kapa0/Kapaq*((1 - Alpha) + Alpha*exp(-Beta*(Kapaq - Kapa0))); - dam = 1. - Kapa0/Kapaq * ((Alpha-Kapaq)/(Alpha - Kapa0)); - dam = std::min(dam,1.); + Real & Equistrain_rate, + Real & Kapaq, + __attribute__((unused)) Real dt, + Real & crit_strain) { + ///Richeton Formula + // Real cst1 = 190e6 - 470e3 * 298.; + // Real cst2 = 2*1.38*1000000*298/5.14; + // Real cst3 = 7.46e15*exp(-1*((90000)/(6.022*1.38*298))); + // Real cst4 = Equistrain_rate/cst3; + // Real cst5 = pow(cst4,1/6.37); + // Real cst6 = log(cst5 + sqrt(1+cst5*cst5)); + // Real cst7 = cst1+cst2*cst6; + // Real cst8 = cst7/(3090.*1e6*1.48); + ///Zhou & Molinari Formula + // Real cst10 = (75./3090.)*(1.+(Equistrain_rate/Beta)); + // Real cst10 = (Kapa0)*(1.+(Equistrain_rate/Beta)); + /// Real cst10 = (75./3090.)*(1.+(Equistrain_rate/(pow(Beta,0.98+1./(1.+Equistrain_rate))))); + // Real cst11 = Alpha; + // Real cst11 = Alpha*(1.+(Equistrain_rate/Beta)); + // Real cst11 = 0; + // if (Equistrain_rate >= 750000.){ + // cst11 = Alpha*(1.+(Equistrain_rate/500000.)); + //} else { + // cst11 = Alpha*(1.+(Equistrain_rate/1000000.)); + //} + + // epsilono + // Real cst10 = Kapa0; + // Real cst10 = (Kapa0)*(1.+(Equistrain_rate/Beta)); + Real cst10 = (75./3090.)*(1.+(Equistrain_rate/Beta)); + // Real cst10 = Equistrain_rate; + + //test + // Real sigmaco=(this->E)*Kapa0*Alpha/(Alpha-Kapa0); + // Real cst10= (sigmaco/this->E)*(1.+(Equistrain_rate/Beta))/(1.+sigmaco/(Alpha*this->E)); + //deg2 + // Real cst10 = (Kapa0)*(1.+(Equistrain_rate/Beta)+(Equistrain_rate/Beta)*(Equistrain_rate/Beta)); + //ln + // Real cst10 = (Kapa0)*(1.+(log(Equistrain_rate+1.)/(Beta/10000.))); + + // epsilonc + Real cst11 = Alpha; + // Real cst11 = Alpha*(1.+(Equistrain_rate/Beta)); + + + + Real cst9 = 0; + if (cst10 > Kapa0) { + if (cst10 < 0.16) { + cst9=cst10; + } else { + cst9=0.16; + } + } else { + cst9 = Kapa0; } - sigma *= 1-dam; + Real Fd1 = Equistrain - cst9; + if (Fd1 > 0) + { + //Real Fd = Equistrain - Kapaq; + // if (Fd > 0) + //{ + Real dam1 = 1. - cst9/Equistrain * ((cst11-Equistrain)/(cst11 - cst9)); + // Real dam1 = 1. - cst9/Equistrain * ((Alpha-Equistrain)/(Alpha - cst9)); + + if (dam1 > dam){ + Kapaq = std::max(Equistrain,Kapaq); + if (dam1 >= 0.9999){ + dam =1.; + } + else { + dam = std::min(dam1,1.); + crit_strain=cst9; + } + } + //} + } + + // Real Alpharates= Alpha/(1.+Equistrain_rate/(6000000.)); new law 2 + // Real Alpharates= 0.6/(1.+Equistrain_rate/(1000000.)); + //new law 3 + // Real Alpharates = 0.6/(1+Equistrain_rate/(pow(1e6,0.98+1./(1+Equistrain_rate)))); + // Real Alpharate=0; + // if (Alpharates <= 0.15){ + // Alpharate=0.15; + // } + // else { + // Alpharate=Alpharates; + // } + // Real dam2 = 1. - cst9/Equistrain * ((Alpharate - Equistrain)/(Alpharate - cst9)); + // if (dam2 >= dam){ + // dam = std::min(dam2,1.); + // } + sigma *= 1-dam; + // crit_strain=cst9; } /* -------------------------------------------------------------------------- */ diff --git a/src/model/solid_mechanics/materials/material_damage/material_vreepeerlings_non_local.hh b/src/model/solid_mechanics/materials/material_damage/material_vreepeerlings_non_local.hh index 1e158dcbd..64390d891 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_vreepeerlings_non_local.hh +++ b/src/model/solid_mechanics/materials/material_damage/material_vreepeerlings_non_local.hh @@ -1,100 +1,106 @@ /** * @file material_vreepeerlings_non_local.hh * @author Cyprien Wolff * @author Nicolas Richart * @date Fri Feb 24 16:01:10 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 "material_vreepeerlings.hh" #include "material_non_local.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_VREEPEERLINGS_NON_LOCAL_HH__ #define __AKANTU_MATERIAL_VREEPEERLINGS_NON_LOCAL_HH__ __BEGIN_AKANTU__ /** * Material VreePeerlings Non local * * parameters in the material files : */ template class WeightFunction = BaseWeightFunction> class MaterialVreePeerlingsNonLocal : public MaterialDamageNonLocal { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef MaterialDamageNonLocal MaterialVreePeerlingsNonLocalParent; MaterialVreePeerlingsNonLocal(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialVreePeerlingsNonLocal() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void initMaterial(); /// constitutive law for all element of a type void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); /// constitutive law virtual void computeNonLocalStress(ElementType el_type, GhostType ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// equivalent strain used to compute the criteria for damage evolution ByElementTypeReal equi_strain; /// non local version of equivalent strain ByElementTypeReal equi_strain_non_local; + /// equivalent strain rate used to compute the criteria for damage evolution + ByElementTypeReal equi_strain_rate; + + /// non local version of equivalent strain rate + ByElementTypeReal equi_strain_rate_non_local; + }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_vreepeerlings_non_local_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_MATERIAL_VREEPEERLINGS_NON_LOCAL_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_damage/material_vreepeerlings_non_local_inline_impl.cc b/src/model/solid_mechanics/materials/material_damage/material_vreepeerlings_non_local_inline_impl.cc index 8a66ca688..738f0dfa5 100644 --- a/src/model/solid_mechanics/materials/material_damage/material_vreepeerlings_non_local_inline_impl.cc +++ b/src/model/solid_mechanics/materials/material_damage/material_vreepeerlings_non_local_inline_impl.cc @@ -1,109 +1,145 @@ /** * @file material_vreepeerlings_non_local_inline_impl.cc * @author Nicolas Richart * @author Cyprien Wolff * @date Fri Jun 15 14:33:40 2012 * * @brief Specialization of the material class for the non-local Vree-Peerlings 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 . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template class WeigthFunction> MaterialVreePeerlingsNonLocal::MaterialVreePeerlingsNonLocal(SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialVreePeerlingsNonLocalParent(model, id), equi_strain("equi-strain", id), - equi_strain_non_local("equi-strain_non_local", id) { + equi_strain_non_local("equi-strain_non_local", id), + equi_strain_rate("equi-strain-rate", id), + equi_strain_rate_non_local("equi-strain-rate_non_local", id) { AKANTU_DEBUG_IN(); this->is_non_local = true; this->initInternalVector(this->equi_strain, 1); this->initInternalVector(this->equi_strain_non_local, 1); + this->initInternalVector(this->equi_strain_rate, 1); + this->initInternalVector(this->equi_strain_rate_non_local, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeigthFunction> void MaterialVreePeerlingsNonLocal::initMaterial() { AKANTU_DEBUG_IN(); this->resizeInternalVector(this->equi_strain); this->resizeInternalVector(this->equi_strain_non_local); + this->resizeInternalVector(this->equi_strain_rate); + this->resizeInternalVector(this->equi_strain_rate_non_local); this->registerNonLocalVariable(this->equi_strain, this->equi_strain_non_local, 1); + this->registerNonLocalVariable(this->equi_strain_rate, this->equi_strain_rate_non_local, 1); MaterialVreePeerlingsNonLocalParent::initMaterial(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeigthFunction> void MaterialVreePeerlingsNonLocal::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real * dam = this->damage(el_type, ghost_type).storage(); Real * equi_straint = equi_strain(el_type, ghost_type).storage(); + Real * equi_straint_rate = equi_strain_rate(el_type, ghost_type).storage(); Real * Kapaq = this->Kapa(el_type, ghost_type).storage(); + Real * crit_strain = this->critical_strain(el_type, ghost_type).storage(); + Real dt = this->model->getTimeStep(); + + Vector & elem_filter = this->element_filter(el_type, ghost_type); + Vector & velocity = this->model->getVelocity(); + Vector & strain_rate_vrplgs = this->strain_rate_vreepeerlings(el_type, ghost_type); + + this->model->getFEM().gradientOnQuadraturePoints(velocity, strain_rate_vrplgs, + spatial_dimension, + el_type, ghost_type, &elem_filter); + + Vector::iterator strain_rate_vrplgs_it = + strain_rate_vrplgs.begin(spatial_dimension, spatial_dimension); + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); + types::Matrix & strain_rate = *strain_rate_vrplgs_it; + MaterialVreePeerlings::computeStressOnQuad(grad_u, sigma, *dam, *equi_straint, - *Kapaq); + *equi_straint_rate, + *Kapaq, + dt, + strain_rate, + *crit_strain); ++dam; ++equi_straint; + ++equi_straint_rate; ++Kapaq; + ++strain_rate_vrplgs_it; + ++crit_strain; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeigthFunction> void MaterialVreePeerlingsNonLocal::computeNonLocalStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); - Real * dam = this->damage(el_type, ghost_type).storage(); - Real * Kapaq = this->Kapa(el_type, ghost_type).storage(); + Real * dam = this->damage(el_type, ghost_type).storage(); + Real * Kapaq = this->Kapa(el_type, ghost_type).storage(); Real * equi_strain_nl = this->equi_strain_non_local(el_type, ghost_type).storage(); + Real * equi_strain_rate_nl = this->equi_strain_rate_non_local(el_type, ghost_type).storage(); + Real dt = this->model->getTimeStep(); + Real * crit_strain = this->critical_strain(el_type, ghost_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); - this->computeDamageAndStressOnQuad(sigma, *dam, *equi_strain_nl, *Kapaq); + this->computeDamageAndStressOnQuad(sigma, *dam, *equi_strain_nl, *equi_strain_rate_nl, *Kapaq, dt, *crit_strain); ++dam; - ++Kapaq; ++equi_strain_nl; + ++equi_strain_rate_nl; + ++Kapaq; + ++crit_strain; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } diff --git a/src/model/solid_mechanics/materials/material_elastic_inline_impl.cc b/src/model/solid_mechanics/materials/material_elastic_inline_impl.cc index 15b6a4c1a..f68fb8cc9 100644 --- a/src/model/solid_mechanics/materials/material_elastic_inline_impl.cc +++ b/src/model/solid_mechanics/materials/material_elastic_inline_impl.cc @@ -1,94 +1,95 @@ /** * @file material_elastic_inline_impl.cc * @author Nicolas Richart * @date Tue Jul 27 11:57:43 2010 * * @brief Implementation of the inline functions of the material 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 . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template inline void MaterialElastic::computeStressOnQuad(const types::Matrix & grad_u, types::Matrix & sigma) { Real trace = grad_u.trace();/// trace = (\nabla u)_{kk} /// \sigma_{ij} = \lambda * (\nabla u)_{kk} * \delta_{ij} + \mu * (\nabla u_{ij} + \nabla u_{ji}) for (UInt i = 0; i < spatial_dimension; ++i) { for (UInt j = 0; j < spatial_dimension; ++j) { sigma(i, j) = (i == j)*lambda*trace + mu*(grad_u(i, j) + grad_u(j, i)); } } + } /* -------------------------------------------------------------------------- */ template<> inline void MaterialElastic<1>::computeStressOnQuad(const types::Matrix & grad_u, types::Matrix & sigma) { sigma(0, 0) = E*grad_u(0, 0); } /* -------------------------------------------------------------------------- */ template void MaterialElastic::computeTangentModuliOnQuad(types::Matrix & tangent) { UInt n = tangent.cols(); //Real Ep = E/((1+nu)*(1-2*nu)); Real Miiii = lambda + 2*mu; Real Miijj = lambda; Real Mijij = mu; if(spatial_dimension == 1) tangent(0, 0) = E; else tangent(0, 0) = Miiii; // test of dimension should by optimized out by the compiler due to the template if(spatial_dimension >= 2) { tangent(1, 1) = Miiii; tangent(0, 1) = Miijj; tangent(1, 0) = Miijj; tangent(n - 1, n - 1) = Mijij; } if(spatial_dimension == 3) { tangent(2, 2) = Miiii; tangent(0, 2) = Miijj; tangent(1, 2) = Miijj; tangent(2, 0) = Miijj; tangent(2, 1) = Miijj; tangent(3, 3) = Mijij; tangent(4, 4) = Mijij; } } /* -------------------------------------------------------------------------- */ template inline Real MaterialElastic::getStableTimeStep(Real h, __attribute__ ((unused)) const Element & element) { return (h/getPushWaveSpeed()); } diff --git a/src/model/solid_mechanics/materials/material_non_local.hh b/src/model/solid_mechanics/materials/material_non_local.hh index 4d54ba221..848387079 100644 --- a/src/model/solid_mechanics/materials/material_non_local.hh +++ b/src/model/solid_mechanics/materials/material_non_local.hh @@ -1,200 +1,215 @@ /** * @file material_non_local.hh * @author Nicolas Richart * @date Thu Jul 28 11:17:41 2011 * * @brief Material class that handle the non locality of a law for example damage. * * @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 "material.hh" #include "aka_grid.hh" #include "fem.hh" #include "weight_function.hh" namespace akantu { class GridSynchronizer; } /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_NON_LOCAL_HH__ #define __AKANTU_MATERIAL_NON_LOCAL_HH__ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template class WeightFunction = BaseWeightFunction> class MaterialNonLocal : public virtual Material { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialNonLocal(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialNonLocal(); template class PairList : public ByElementType > {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// read properties virtual bool parseParam(const std::string & key, const std::string & value, const ID & id); /// initialize the material computed parameter virtual void initMaterial(); // void computeQuadraturePointsNeighborhoudVolumes(ByElementTypeReal & volumes) const; virtual void updateResidual(GhostType ghost_type); virtual void computeAllNonLocalStresses(GhostType ghost_type = _not_ghost); // void removeDamaged(const ByElementTypeReal & damage, Real thresold); void savePairs(const std::string & filename) const; void neighbourhoodStatistics(const std::string & filename) const; protected: void updatePairList(const ByElementTypeReal & quadrature_points_coordinates); void computeWeights(const ByElementTypeReal & quadrature_points_coordinates); void createCellList(ByElementTypeReal & quadrature_points_coordinates); void fillCellList(const ByElementTypeReal & quadrature_points_coordinates, const GhostType & ghost_type); /// constitutive law - virtual void computeNonLocalStress(GhostType ghost_type = _not_ghost) = 0; + virtual void computeNonLocalStresses(GhostType ghost_type = _not_ghost) = 0; template void weightedAvergageOnNeighbours(const ByElementTypeVector & to_accumulate, ByElementTypeVector & accumulated, UInt nb_degree_of_freedom, GhostType ghost_type2 = _not_ghost) const; // template // void accumulateOnNeighbours(const ByElementTypeVector & to_accumulate, // ByElementTypeVector & accumulated, // UInt nb_degree_of_freedom) const; virtual inline UInt getNbDataToPack(const Element & element, SynchronizationTag tag) const; virtual inline UInt getNbDataToUnpack(const Element & element, SynchronizationTag tag) const; virtual inline void packData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag) const; virtual inline void unpackData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag); + virtual UInt getNbDataToPack(SynchronizationTag tag) const { return Material::getNbDataToPack(tag); } + virtual UInt getNbDataToUnpack(SynchronizationTag tag) const { return Material::getNbDataToUnpack(tag); } + + virtual void packData(CommunicationBuffer & buffer, + const UInt index, + SynchronizationTag tag) const { + Material::packData(buffer, index, tag); + } + + virtual void unpackData(CommunicationBuffer & buffer, + const UInt index, + SynchronizationTag tag) { + Material::unpackData(buffer, index, tag); + } + virtual inline void onElementsAdded(const Vector & element_list); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: void registerNonLocalVariable(ByElementTypeReal & local, ByElementTypeReal & non_local, UInt nb_degree_of_freedom) { ID id = local.getID(); NonLocalVariable & non_local_variable = non_local_variables[id]; non_local_variable.local_variable = &local; non_local_variable.non_local_variable = &non_local; non_local_variable.non_local_variable_nb_component = nb_degree_of_freedom; } AKANTU_GET_MACRO(PairList, pair_list, const PairList &) AKANTU_GET_MACRO(Radius, radius, Real); AKANTU_GET_MACRO(CellList, *cell_list, const RegularGrid &) /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// the non local radius Real radius; /// the weight function used WeightFunction * weight_func; private: /// the pairs of quadrature points PairList pair_list; /// the weights associated to the pairs PairList pair_weight; /// the regular grid to construct/update the pair lists RegularGrid * cell_list; /// the types of the existing pairs typedef std::set< std::pair > pair_type; pair_type existing_pairs[2]; /// specify if the weights should be updated and at which rate UInt update_weights; /// count the number of calls of computeStress UInt compute_stress_calls; struct NonLocalVariable { ByElementTypeVector * local_variable; ByElementTypeVector * non_local_variable; UInt non_local_variable_nb_component; }; std::map non_local_variables; bool is_creating_grid; GridSynchronizer * grid_synchronizer; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_non_local_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_MATERIAL_NON_LOCAL_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_non_local_inline_impl.cc b/src/model/solid_mechanics/materials/material_non_local_inline_impl.cc index 13b1e1a6b..99e7af095 100644 --- a/src/model/solid_mechanics/materials/material_non_local_inline_impl.cc +++ b/src/model/solid_mechanics/materials/material_non_local_inline_impl.cc @@ -1,761 +1,761 @@ /** * @file material_non_local_inline_impl.cc * @author Nicolas Richart * @date Thu Aug 25 11:59:39 2011 * * @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 . * */ /* -------------------------------------------------------------------------- */ __END_AKANTU__ /* -------------------------------------------------------------------------- */ #include "aka_types.hh" #include "grid_synchronizer.hh" #include "synchronizer_registry.hh" #include "integrator.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template class WeightFunction> MaterialNonLocal::MaterialNonLocal(SolidMechanicsModel & model, const ID & id) : Material(model, id), radius(100.), weight_func(NULL), cell_list(NULL), update_weights(0), compute_stress_calls(0), is_creating_grid(false), grid_synchronizer(NULL) { AKANTU_DEBUG_IN(); this->registerParam("radius" , radius , 100., _pat_readable , "Non local radius"); this->registerParam("UpdateWeights" , update_weights, 0U , _pat_modifiable, "Update weights frequency"); this->registerParam("Weight function", weight_func , _pat_internal); this->is_non_local = true; this->weight_func = new WeightFunction(*this); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> MaterialNonLocal::~MaterialNonLocal() { AKANTU_DEBUG_IN(); delete cell_list; delete weight_func; delete grid_synchronizer; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::initMaterial() { AKANTU_DEBUG_IN(); // Material::initMaterial(); Mesh & mesh = this->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); createCellList(quadrature_points_coordinates); updatePairList(quadrature_points_coordinates); computeWeights(quadrature_points_coordinates); weight_func->setRadius(radius); weight_func->init(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::updateResidual(GhostType ghost_type) { AKANTU_DEBUG_IN(); // Update the weights for the non local variable averaging if(ghost_type == _not_ghost && this->update_weights && (this->compute_stress_calls % this->update_weights == 0)) { ByElementTypeReal quadrature_points_coordinates("quadrature_points_coordinates", id); Mesh & mesh = this->model->getFEM().getMesh(); mesh.initByElementTypeVector(quadrature_points_coordinates, spatial_dimension, 0); computeQuadraturePointsCoordinates(quadrature_points_coordinates, _not_ghost); computeQuadraturePointsCoordinates(quadrature_points_coordinates, _ghost); computeWeights(quadrature_points_coordinates); } if(ghost_type == _not_ghost) ++this->compute_stress_calls; computeAllStresses(ghost_type); - computeNonLocalStress(ghost_type); + computeNonLocalStresses(ghost_type); assembleResidual(ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::computeAllNonLocalStresses(GhostType ghost_type) { // Update the weights for the non local variable averaging if(ghost_type == _not_ghost) { if(this->update_weights && (this->compute_stress_calls % this->update_weights == 0)) { this->model->getSynchronizerRegistry().asynchronousSynchronize(_gst_mnl_weight); ByElementTypeReal quadrature_points_coordinates("quadrature_points_coordinates", id); Mesh & mesh = this->model->getFEM().getMesh(); mesh.initByElementTypeVector(quadrature_points_coordinates, spatial_dimension, 0); computeQuadraturePointsCoordinates(quadrature_points_coordinates, _not_ghost); computeQuadraturePointsCoordinates(quadrature_points_coordinates, _ghost); this->model->getSynchronizerRegistry().waitEndSynchronize(_gst_mnl_weight); computeWeights(quadrature_points_coordinates); } typename std::map::iterator it = non_local_variables.begin(); typename std::map::iterator end = non_local_variables.end(); for(;it != end; ++it) { NonLocalVariable & non_local_variable = it->second; resizeInternalVector(*non_local_variable.non_local_variable); this->weightedAvergageOnNeighbours(*non_local_variable.local_variable, *non_local_variable.non_local_variable, non_local_variable.non_local_variable_nb_component, _not_ghost); } ++this->compute_stress_calls; } else { typename std::map::iterator it = non_local_variables.begin(); typename std::map::iterator end = non_local_variables.end(); for(;it != end; ++it) { NonLocalVariable & non_local_variable = it->second; this->weightedAvergageOnNeighbours(*non_local_variable.local_variable, *non_local_variable.non_local_variable, non_local_variable.non_local_variable_nb_component, _ghost); } - computeNonLocalStress(_not_ghost); + computeNonLocalStresses(_not_ghost); } } /* -------------------------------------------------------------------------- */ template class WeightFunction> template void MaterialNonLocal::weightedAvergageOnNeighbours(const ByElementTypeVector & to_accumulate, ByElementTypeVector & accumulated, UInt nb_degree_of_freedom, GhostType ghost_type2) const { AKANTU_DEBUG_IN(); UInt existing_pairs_num = 0; if (ghost_type2 == _ghost) existing_pairs_num = 1; pair_type::const_iterator first_pair_types = existing_pairs[existing_pairs_num].begin(); pair_type::const_iterator last_pair_types = existing_pairs[existing_pairs_num].end(); GhostType ghost_type1 = _not_ghost; // does not make sens the ghost vs ghost so this should always by _not_ghost for (; first_pair_types != last_pair_types; ++first_pair_types) { const Vector & pairs = pair_list(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); const Vector & weights = pair_weight(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); const Vector & to_acc = to_accumulate(first_pair_types->second, ghost_type2); Vector & acc = accumulated(first_pair_types->first, ghost_type1); Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); Vector::const_iterator< types::Vector > pair_w = weights.begin(2); for(;first_pair != last_pair; ++first_pair, ++pair_w) { UInt q1 = (*first_pair)(0); UInt q2 = (*first_pair)(1); for(UInt d = 0; d < nb_degree_of_freedom; ++d){ acc(q1, d) += (*pair_w)(0) * to_acc(q2, d); if(ghost_type2 != _ghost) acc(q2, d) += (*pair_w)(1) * to_acc(q1, d); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::createCellList(ByElementTypeReal & quadrature_points_coordinates) { AKANTU_DEBUG_IN(); const Real safety_factor = 1.2; // for the cell grid spacing Mesh & mesh = this->model->getFEM().getMesh(); mesh.computeBoundingBox(); Real lower_bounds[spatial_dimension]; Real upper_bounds[spatial_dimension]; mesh.getLocalLowerBounds(lower_bounds); mesh.getLocalUpperBounds(upper_bounds); Real spacing[spatial_dimension]; for (UInt i = 0; i < spatial_dimension; ++i) { spacing[i] = radius * safety_factor; } cell_list = new RegularGrid(spatial_dimension, lower_bounds, upper_bounds, spacing); fillCellList(quadrature_points_coordinates, _not_ghost); is_creating_grid = true; SynchronizerRegistry & synch_registry = this->model->getSynchronizerRegistry(); std::stringstream sstr; sstr << id << ":grid_synchronizer"; grid_synchronizer = GridSynchronizer::createGridSynchronizer(mesh, *cell_list, sstr.str()); synch_registry.registerSynchronizer(*grid_synchronizer, _gst_mnl_for_average); synch_registry.registerSynchronizer(*grid_synchronizer, _gst_mnl_weight); is_creating_grid = false; this->computeQuadraturePointsCoordinates(quadrature_points_coordinates, _ghost); fillCellList(quadrature_points_coordinates, _ghost); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::fillCellList(const ByElementTypeReal & quadrature_points_coordinates, const GhostType & ghost_type) { Mesh & mesh = this->model->getFEM().getMesh(); QuadraturePoint q; Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type); q.ghost_type = ghost_type; for(; it != last_type; ++it) { Vector & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_quad = this->model->getFEM().getNbQuadraturePoints(*it, ghost_type); const Vector & quads = quadrature_points_coordinates(*it, ghost_type); q.type = *it; Vector::const_iterator quad = quads.begin(spatial_dimension); UInt * elem = elem_filter.storage(); for (UInt e = 0; e < nb_element; ++e) { q.element = *elem; for (UInt nq = 0; nq < nb_quad; ++nq) { q.num_point = nq; q.setPosition(*quad); cell_list->insert(q, *quad); ++quad; } ++elem; } } } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::updatePairList(const ByElementTypeReal & quadrature_points_coordinates) { AKANTU_DEBUG_IN(); Mesh & mesh = this->model->getFEM().getMesh(); GhostType ghost_type = _not_ghost; // generate the pair of neighbor depending of the cell_list 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) { // Preparing datas const Vector & quads = quadrature_points_coordinates(*it, ghost_type); Vector::const_iterator first_quad = quads.begin(spatial_dimension); Vector::const_iterator last_quad = quads.end(spatial_dimension); ByElementTypeUInt & pairs = pair_list(ByElementTypeUInt("pairs", id, memory_id), *it, ghost_type); ElementType current_element_type = _not_defined; GhostType current_ghost_type = _casper; UInt existing_pairs_num = 0; Vector * neighbors = NULL; UInt my_num_quad = 0; // loop over quad points for(;first_quad != last_quad; ++first_quad, ++my_num_quad) { RegularGrid::Cell cell = cell_list->getCell(*first_quad); RegularGrid::neighbor_cells_iterator first_neigh_cell = cell_list->beginNeighborCells(cell); RegularGrid::neighbor_cells_iterator last_neigh_cell = cell_list->endNeighborCells(cell); // loop over neighbors cells of the one containing the current quadrature // point for (; first_neigh_cell != last_neigh_cell; ++first_neigh_cell) { RegularGrid::iterator first_neigh_quad = cell_list->beginCell(*first_neigh_cell); RegularGrid::iterator last_neigh_quad = cell_list->endCell(*first_neigh_cell); // loop over the quadrature point in the current cell of the cell list for (;first_neigh_quad != last_neigh_quad; ++first_neigh_quad){ QuadraturePoint & quad = *first_neigh_quad; UInt nb_quad_per_elem = this->model->getFEM().getNbQuadraturePoints(quad.type, quad.ghost_type); UInt neigh_num_quad = quad.element * nb_quad_per_elem + quad.num_point; // little optimization to not search in the map at each quad points if(quad.type != current_element_type || quad.ghost_type != current_ghost_type) { // neigh_quad_positions = quadrature_points_coordinates(quad.type, // quad.ghost_type).storage(); current_element_type = quad.type; current_ghost_type = quad.ghost_type; existing_pairs_num = quad.ghost_type == _not_ghost ? 0 : 1; if(!pairs.exists(current_element_type, current_ghost_type)) { neighbors = &(pairs.alloc(0, 2, current_element_type, current_ghost_type)); } else { neighbors = &(pairs(current_element_type, current_ghost_type)); } existing_pairs[existing_pairs_num].insert(std::pair(*it, current_element_type)); } // types::RVector neigh_quad(neigh_quad_positions + neigh_num_quad * spatial_dimension, // spatial_dimension); const types::RVector & neigh_quad = quad.getPosition(); Real distance = first_quad->distance(neigh_quad); if(distance <= radius && (current_ghost_type == _ghost || (current_ghost_type == _not_ghost && my_num_quad <= neigh_num_quad))) { // sotring only half lists UInt pair[2]; pair[0] = my_num_quad; pair[1] = neigh_num_quad; neighbors->push_back(pair); } // } } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::computeWeights(const ByElementTypeReal & quadrature_points_coordinates) { AKANTU_DEBUG_IN(); GhostType ghost_type1; ghost_type1 = _not_ghost; ByElementTypeReal quadrature_points_volumes("quadrature_points_volumes", id, memory_id); this->model->getFEM().getMesh().initByElementTypeVector(quadrature_points_volumes, 1, 0); resizeInternalVector(quadrature_points_volumes); const ByElementTypeReal & jacobians_by_type = this->model->getFEM().getIntegratorInterface().getJacobians(); weight_func->updateInternals(quadrature_points_volumes); for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType) gt; UInt existing_pairs_num = gt - _not_ghost; pair_type::iterator first_pair_types = existing_pairs[existing_pairs_num].begin(); pair_type::iterator last_pair_types = existing_pairs[existing_pairs_num].end(); // Compute the weights for (; first_pair_types != last_pair_types; ++first_pair_types) { ElementType type1 = first_pair_types->first; ElementType type2 = first_pair_types->second; const Vector & pairs = pair_list(type1, ghost_type1)(type2, ghost_type2); std::string ghost_id = ""; if (ghost_type1 == _ghost) ghost_id = ":ghost"; ByElementTypeReal & weights_type_1 = pair_weight(type1, ghost_type1); std::stringstream sstr; sstr << id << ":pair_weight:" << type1 << ghost_id; weights_type_1.setID(sstr.str()); Vector * tmp_weight = NULL; if(!weights_type_1.exists(type2, ghost_type2)) { tmp_weight = &(weights_type_1.alloc(0, 2, type2, ghost_type2)); } else { tmp_weight = &(weights_type_1(type2, ghost_type2)); } Vector & weights = *tmp_weight; weights.resize(pairs.getSize()); weights.clear(); const Vector & jacobians_1 = jacobians_by_type(type1, ghost_type1); const Vector & jacobians_2 = jacobians_by_type(type2, ghost_type2); // const Vector & elem_filter1 = element_filter(type1, ghost_type1); // const Vector & elem_filter2 = element_filter(type2, ghost_type2); UInt nb_quad1 = this->model->getFEM().getNbQuadraturePoints(type1); UInt nb_quad2 = this->model->getFEM().getNbQuadraturePoints(type2); // UInt nb_tot_quad1 = nb_quad1 * elem_filter1.getSize();; // UInt nb_tot_quad2 = nb_quad2 * elem_filter2.getSize();; Vector & quads_volumes1 = quadrature_points_volumes(type1, ghost_type1); Vector & quads_volumes2 = quadrature_points_volumes(type2, ghost_type2); Vector::const_iterator iquads1; Vector::const_iterator iquads2; iquads1 = quadrature_points_coordinates(type1, ghost_type1).begin(spatial_dimension); iquads2 = quadrature_points_coordinates(type2, ghost_type2).begin(spatial_dimension); Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); Vector::iterator weight = weights.begin(2); this->weight_func->selectType(type1, ghost_type1, type2, ghost_type2); // Weight function for(;first_pair != last_pair; ++first_pair, ++weight) { UInt _q1 = (*first_pair)(0); UInt _q2 = (*first_pair)(1); const types::RVector & pos1 = iquads1[_q1]; const types::RVector & pos2 = iquads2[_q2]; QuadraturePoint q1(_q1 / nb_quad1, _q1 % nb_quad1, _q1, pos1, type1, ghost_type1); QuadraturePoint q2(_q2 / nb_quad2, _q2 % nb_quad2, _q2, pos2, type2, ghost_type2); Real r = pos1.distance(pos2); Real w2J2 = jacobians_2(_q2); (*weight)(0) = w2J2 * this->weight_func->operator()(r, q1, q2); if(ghost_type2 != _ghost && _q1 != _q2) { Real w1J1 = jacobians_1(_q1); (*weight)(1) = w1J1 * this->weight_func->operator()(r, q2, q1); } else (*weight)(1) = 0; quads_volumes1(_q1) += (*weight)(0); if(ghost_type2 != _ghost) quads_volumes2(_q2) += (*weight)(1); } } } //normalize the weights for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType) gt; UInt existing_pairs_num = gt - _not_ghost; pair_type::iterator first_pair_types = existing_pairs[existing_pairs_num].begin(); pair_type::iterator last_pair_types = existing_pairs[existing_pairs_num].end(); first_pair_types = existing_pairs[existing_pairs_num].begin(); for (; first_pair_types != last_pair_types; ++first_pair_types) { ElementType type1 = first_pair_types->first; ElementType type2 = first_pair_types->second; const Vector & pairs = pair_list(type1, ghost_type1)(type2, ghost_type2); Vector & weights = pair_weight(type1, ghost_type1)(type2, ghost_type2); Vector & quads_volumes1 = quadrature_points_volumes(type1, ghost_type1); Vector & quads_volumes2 = quadrature_points_volumes(type2, ghost_type2); Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); Vector::iterator weight = weights.begin(2); for(;first_pair != last_pair; ++first_pair, ++weight) { UInt q1 = (*first_pair)(0); UInt q2 = (*first_pair)(1); (*weight)(0) *= 1. / quads_volumes1(q1); if(ghost_type2 != _ghost) (*weight)(1) *= 1. / quads_volumes2(q2); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> bool MaterialNonLocal::parseParam(const std::string & key, const std::string & value, __attribute__((unused)) const ID & id) { std::stringstream sstr(value); if(key == "radius") { sstr >> radius; } else if(key == "UpdateWeights") { sstr >> update_weights; } else if(!weight_func->parseParam(key, value)) return Material::parseParam(key, value, id); return true; } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::savePairs(const std::string & filename) const { std::ofstream pout; std::stringstream sstr; StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Int prank = comm.whoAmI(); sstr << filename << "." << prank; pout.open(sstr.str().c_str()); GhostType ghost_type1; ghost_type1 = _not_ghost; for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType) gt; UInt existing_pairs_num = gt - _not_ghost; pair_type::const_iterator first_pair_types = existing_pairs[existing_pairs_num].begin(); pair_type::const_iterator last_pair_types = existing_pairs[existing_pairs_num].end(); for (; first_pair_types != last_pair_types; ++first_pair_types) { const Vector & pairs = pair_list(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); const Vector & weights = pair_weight(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); pout << "Types : " << first_pair_types->first << " (" << ghost_type1 << ") - " << first_pair_types->second << " (" << ghost_type2 << ")" << std::endl; Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); Vector::const_iterator pair_w = weights.begin(2); for(;first_pair != last_pair; ++first_pair, ++pair_w) { UInt q1 = (*first_pair)(0); UInt q2 = (*first_pair)(1); pout << q1 << " " << q2 << " "<< (*pair_w)(0) << " " << (*pair_w)(1) << std::endl; } } } } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::neighbourhoodStatistics(const std::string & filename) const { std::ofstream pout; pout.open(filename.c_str()); const Mesh & mesh = this->model->getFEM().getMesh(); GhostType ghost_type1; ghost_type1 = _not_ghost; StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Int prank = comm.whoAmI(); for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType) gt; UInt existing_pairs_num = gt - _not_ghost; ByElementTypeUInt nb_neighbors("nb_neighbours", id, memory_id); mesh.initByElementTypeVector(nb_neighbors, 1, spatial_dimension); resizeInternalVector(nb_neighbors); pair_type::const_iterator first_pair_types = existing_pairs[existing_pairs_num].begin(); pair_type::const_iterator last_pair_types = existing_pairs[existing_pairs_num].end(); for (; first_pair_types != last_pair_types; ++first_pair_types) { const Vector & pairs = pair_list(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); if(prank == 0) { pout << ghost_type2 << " "; pout << "Types : " << first_pair_types->first << " " << first_pair_types->second << std::endl; } Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); Vector & nb_neigh_1 = nb_neighbors(first_pair_types->first, ghost_type1); Vector & nb_neigh_2 = nb_neighbors(first_pair_types->second, ghost_type2); for(;first_pair != last_pair; ++first_pair) { UInt q1 = (*first_pair)(0); UInt q2 = (*first_pair)(1); ++(nb_neigh_1(q1)); if(q1 != q2) ++(nb_neigh_2(q2)); } } Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type1); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type1); UInt nb_quads = 0; Real sum_nb_neig = 0; UInt max_nb_neig = 0; UInt min_nb_neig = std::numeric_limits::max(); for(; it != last_type; ++it) { Vector & nb_neighor = nb_neighbors(*it, ghost_type1); Vector::iterator nb_neigh = nb_neighor.begin(); Vector::iterator end_neigh = nb_neighor.end(); for (; nb_neigh != end_neigh; ++nb_neigh, ++nb_quads) { UInt nb = *nb_neigh; sum_nb_neig += nb; max_nb_neig = std::max(max_nb_neig, nb); min_nb_neig = std::min(min_nb_neig, nb); } } comm.allReduce(&nb_quads, 1, _so_sum); comm.allReduce(&sum_nb_neig, 1, _so_sum); comm.allReduce(&max_nb_neig, 1, _so_max); comm.allReduce(&min_nb_neig, 1, _so_min); if(prank == 0) { pout << ghost_type2 << " "; pout << "Nb quadrature points: " << nb_quads << std::endl; Real mean_nb_neig = sum_nb_neig / Real(nb_quads); pout << ghost_type2 << " "; pout << "Average nb neighbors: " << mean_nb_neig << "(" << sum_nb_neig << ")" << std::endl; pout << ghost_type2 << " "; pout << "Max nb neighbors: " << max_nb_neig << std::endl; pout << ghost_type2 << " "; pout << "Min nb neighbors: " << min_nb_neig << std::endl; } } pout.close(); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template class WeightFunction> inline UInt MaterialNonLocal::getNbDataToPack(const Element & element, SynchronizationTag tag) const { UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(element.type); if(tag == _gst_mnl_for_average) { typename std::map::const_iterator it = non_local_variables.begin(); typename std::map::const_iterator end = non_local_variables.end(); UInt size = 0; for(;it != end; ++it) { const NonLocalVariable & non_local_variable = it->second; size += non_local_variable.non_local_variable_nb_component; } return size * sizeof(Real) * nb_quadrature_points; } else if(tag == _gst_mnl_weight) return weight_func->getNbData(element, tag); return 0; } /* -------------------------------------------------------------------------- */ template class WeightFunction> inline UInt MaterialNonLocal::getNbDataToUnpack(const Element & element, SynchronizationTag tag) const { UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(element.type); if(tag == _gst_mnl_for_average) { typename std::map::const_iterator it = non_local_variables.begin(); typename std::map::const_iterator end = non_local_variables.end(); UInt size = 0; for(;it != end; ++it) { const NonLocalVariable & non_local_variable = it->second; size += non_local_variable.non_local_variable_nb_component; } return size * sizeof(Real) * nb_quadrature_points; } else if(tag == _gst_mnl_weight) return weight_func->getNbData(element, tag); return 0; } /* -------------------------------------------------------------------------- */ template class WeightFunction> inline void MaterialNonLocal::packData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag) const { if(tag == _gst_mnl_for_average) { UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(element.type); typename std::map::const_iterator it = non_local_variables.begin(); typename std::map::const_iterator end = non_local_variables.end(); for(;it != end; ++it) { const NonLocalVariable & non_local_variable = it->second; Vector::iterator local_var = (*non_local_variable.local_variable)(element.type, _not_ghost).begin(non_local_variable.non_local_variable_nb_component); local_var += element.element * nb_quadrature_points; for (UInt q = 0; q < nb_quadrature_points; ++q, ++local_var) buffer << *local_var; } } else if(tag == _gst_mnl_weight) return weight_func->packData(buffer, element, tag); } /* -------------------------------------------------------------------------- */ template class WeightFunction> inline void MaterialNonLocal::unpackData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag) { if(tag == _gst_mnl_for_average) { UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(element.type); typename std::map::iterator it = non_local_variables.begin(); typename std::map::iterator end = non_local_variables.end(); for(;it != end; ++it) { NonLocalVariable & non_local_variable = it->second; Vector::iterator local_var = (*non_local_variable.local_variable)(element.type, _ghost).begin(non_local_variable.non_local_variable_nb_component); local_var += element.element * nb_quadrature_points; for (UInt q = 0; q < nb_quadrature_points; ++q, ++local_var) buffer >> *local_var; } } else if(tag == _gst_mnl_weight) return weight_func->unpackData(buffer, element, tag); } /* -------------------------------------------------------------------------- */ template class WeightFunction> inline void MaterialNonLocal::onElementsAdded(const Vector & element_list) { AKANTU_DEBUG_IN(); if(is_creating_grid) { Int my_index = model->getInternalIndexFromID(id); std::cout << "Toto :" << my_index << std::endl; AKANTU_DEBUG_ASSERT(my_index != -1, "Something horrible happen, the model does not know the material " << id); Vector::const_iterator it = element_list.begin(); Vector::const_iterator end = element_list.end(); for (; it != end; ++it) { const Element & elem = *it; model->getElementIndexByMaterial(elem.type, elem.ghost_type)(elem.element) = this->addElement(elem.type, elem.element, elem.ghost_type); model->getElementMaterial(elem.type, elem.ghost_type)(elem.element) = UInt(my_index); } } Material::onElementsAdded(element_list); AKANTU_DEBUG_OUT(); } diff --git a/src/model/solid_mechanics/materials/weight_function.hh b/src/model/solid_mechanics/materials/weight_function.hh index 0d20a4c40..0e10aec12 100644 --- a/src/model/solid_mechanics/materials/weight_function.hh +++ b/src/model/solid_mechanics/materials/weight_function.hh @@ -1,335 +1,415 @@ /** * @file weight_function.hh * @author Nicolas Richart * @date Thu Mar 8 16:17:00 2012 * * @brief Weight functions for non local materials * * @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_types.hh" #include "solid_mechanics_model.hh" +#include /* -------------------------------------------------------------------------- */ #include #ifndef __AKANTU_WEIGHT_FUNCTION_HH__ #define __AKANTU_WEIGHT_FUNCTION_HH__ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* Normal weight function */ /* -------------------------------------------------------------------------- */ template class BaseWeightFunction { public: BaseWeightFunction(Material & material) : material(material) {} virtual ~BaseWeightFunction() {} virtual void init() {}; virtual void updateInternals(__attribute__((unused)) const ByElementTypeReal & quadrature_points_coordinates) {}; /* ------------------------------------------------------------------------ */ inline void setRadius(Real radius) { R = radius; R2 = R * R; } /* ------------------------------------------------------------------------ */ inline void selectType(__attribute__((unused)) ElementType type1, __attribute__((unused)) GhostType ghost_type1, __attribute__((unused)) ElementType type2, __attribute__((unused)) GhostType ghost_type2) { } /* ------------------------------------------------------------------------ */ inline Real operator()(Real r, __attribute__((unused)) QuadraturePoint & q1, __attribute__((unused)) QuadraturePoint & q2) { Real w = 0; if(r <= R) { Real alpha = (1. - r*r / R2); w = alpha * alpha; // *weight = 1 - sqrt(r / radius); } return w; } /* ------------------------------------------------------------------------ */ bool parseParam(__attribute__((unused)) const std::string & key, __attribute__((unused)) const std::string & value) { return false; } /* ------------------------------------------------------------------------ */ virtual void printself(std::ostream & stream) const { stream << "BaseWeightFunction"; } public: virtual UInt getNbData(__attribute__((unused)) const Element & element, __attribute__((unused)) SynchronizationTag tag) const { return 0; } - virtual inline void packData(CommunicationBuffer & buffer, - const Element & element, - SynchronizationTag tag) const {} + virtual inline void packData(__attribute__((unused)) CommunicationBuffer & buffer, + __attribute__((unused)) const Element & element, + __attribute__((unused)) SynchronizationTag tag) const {} - virtual inline void unpackData(CommunicationBuffer & buffer, - const Element & element, - SynchronizationTag tag) {} + virtual inline void unpackData(__attribute__((unused)) CommunicationBuffer & buffer, + __attribute__((unused)) const Element & element, + __attribute__((unused)) SynchronizationTag tag) {} protected: Material & material; Real R; Real R2; }; /* -------------------------------------------------------------------------- */ /* Damage weight function */ /* -------------------------------------------------------------------------- */ template class DamagedWeightFunction : public BaseWeightFunction { public: DamagedWeightFunction(Material & material) : BaseWeightFunction(material) {} inline void selectType(__attribute__((unused)) ElementType type1, __attribute__((unused)) GhostType ghost_type1, ElementType type2, GhostType ghost_type2) { selected_damage = &this->material.getVector("damage", type2, ghost_type2); } inline Real operator()(Real r, __attribute__((unused)) QuadraturePoint & q1, QuadraturePoint & q2) { UInt quad = q2.global_num; Real D = (*selected_damage)(quad); Real Radius = (1.-D)*(1.-D) * this->R2; if(Radius < Math::getTolerance()) { Radius = 0.01 * 0.01 * this->R2; } Real alpha = std::max(0., 1. - r*r / Radius); Real w = alpha * alpha; return w; } /* ------------------------------------------------------------------------ */ bool parseParam(__attribute__((unused)) const std::string & key, __attribute__((unused)) const std::string & value) { return false; } /* ------------------------------------------------------------------------ */ virtual void printself(std::ostream & stream) const { stream << "DamagedWeightFunction"; } private: const Vector * selected_damage; }; /* -------------------------------------------------------------------------- */ /* Remove damaged weight function */ /* -------------------------------------------------------------------------- */ template class RemoveDamagedWeightFunction : public BaseWeightFunction { public: RemoveDamagedWeightFunction(Material & material) : BaseWeightFunction(material) {} inline void selectType(__attribute__((unused)) ElementType type1, __attribute__((unused)) GhostType ghost_type1, ElementType type2, GhostType ghost_type2) { selected_damage = &this->material.getVector("damage", type2, ghost_type2); } inline Real operator()(Real r, __attribute__((unused)) QuadraturePoint & q1, QuadraturePoint & q2) { UInt quad = q2.global_num; if(q1.global_num == quad) return 1.; Real D = (*selected_damage)(quad); Real w = 0.; if(D < damage_limit) { - Real alpha = std::max(0., 1. - r*r / this->R2); - w = alpha * alpha; + + ////alpha2beta6//// + // Real alpha = std::max(0., 1. - r*r / this->R2); + // w = alpha * alpha * alpha * alpha * alpha * alpha; + + ////alpha4beta6//// + //Real alpha = std::max(0., 1. - (r*r / this->R2)*(r*r / this->R2)); + //w = alpha * alpha * alpha * alpha * alpha * alpha; + + ////alpha6beta6//// + //Real alpha = std::max(0., 1. - (r*r / this->R2)*(r*r / this->R2)*(r*r / this->R2)); + //w = alpha * alpha * alpha * alpha * alpha * alpha; + + ////alpha8beta6//// + //Real alpha = std::max(0., 1. - (r*r / this->R2)*(r*r / this->R2)*(r*r / this->R2)*(r*r / this->R2)); + //w = alpha * alpha * alpha * alpha * alpha * alpha; + + ////alpha10beta6//// + //Real alpha = std::max(0., 1. - (r*r / this->R2)*(r*r / this->R2)*(r*r / this->R2)*(r*r / this->R2)*(r*r / this->R2)); + //w = alpha * alpha * alpha * alpha * alpha * alpha; + + ////alpha2beta2//// + Real alpha = std::max(0., 1. - r*r / this->R2); + w = alpha * alpha; } return w; } /* ------------------------------------------------------------------------ */ bool parseParam(const std::string & key, const std::string & value) { std::stringstream sstr(value); if(key == "damage_limit") { sstr >> damage_limit; } else return BaseWeightFunction::parseParam(key, value); return true; } /* ------------------------------------------------------------------------ */ virtual void printself(std::ostream & stream) const { stream << "RemoveDamagedWeightFunction [damage_limit: " << damage_limit << "]"; } virtual UInt getNbData(const Element & element, __attribute__((unused)) SynchronizationTag tag) const { UInt nb_quadrature_points = this->material.getModel().getFEM().getNbQuadraturePoints(element.type); UInt size = sizeof(Real) * nb_quadrature_points; return size; } virtual inline void packData(CommunicationBuffer & buffer, const Element & element, __attribute__((unused)) SynchronizationTag tag) const { UInt nb_quadrature_points = this->material.getModel().getFEM().getNbQuadraturePoints(element.type); const Vector & dam_vect = this->material.getVector("damage", element.type, _not_ghost); Vector::const_iterator damage = dam_vect.begin(); damage += element.element * nb_quadrature_points; for (UInt q = 0; q < nb_quadrature_points; ++q, ++damage) buffer << *damage; } virtual inline void unpackData(CommunicationBuffer & buffer, const Element & element, __attribute__((unused)) SynchronizationTag tag) { UInt nb_quadrature_points = this->material.getModel().getFEM().getNbQuadraturePoints(element.type); Vector::iterator damage = this->material.getVector("damage", element.type, _ghost).begin(); damage += element.element * nb_quadrature_points; for (UInt q = 0; q < nb_quadrature_points; ++q, ++damage) buffer >> *damage; } private: /// limit at which a point is considered as complitely broken Real damage_limit; /// internal pointer to the current damage vector const Vector * selected_damage; }; +/* -------------------------------------------------------------------------- */ +/* Remove damaged with damage rate weight function */ +/* -------------------------------------------------------------------------- */ + +template +class RemoveDamagedWithDamageRateWeightFunction : public BaseWeightFunction { +public: + RemoveDamagedWithDamageRateWeightFunction(Material & material) : BaseWeightFunction(material) {} + + inline void selectType(__attribute__((unused)) ElementType type1, + __attribute__((unused)) GhostType ghost_type1, + ElementType type2, + GhostType ghost_type2) { + selected_damage_with_damage_rate = &(this->material.getVector("damage",type2, ghost_type2)); + selected_damage_rate_with_damage_rate = &(this->material.getVector("damage-rate",type2, ghost_type2)); + } + + inline Real operator()(Real r, __attribute__((unused)) QuadraturePoint & q1, QuadraturePoint & q2) { + UInt quad = q2.global_num; + + if(q1.global_num == quad) return 1.; + + Real D = (*selected_damage_with_damage_rate)(quad); + Real w = 0.; + Real alphaexp = 1.; + Real betaexp = 2.; + if(D < damage_limit_with_damage_rate) { + Real alpha = std::max(0., 1. - pow((r*r / this->R2),alphaexp)); + w = pow(alpha, betaexp); + } + + return w; + } + /* ------------------------------------------------------------------------ */ + bool setParam(const std::string & key, + const std::string & value) { + std::stringstream sstr(value); + if(key == "damage_limit") { sstr >> damage_limit_with_damage_rate; } + else return BaseWeightFunction::setParam(key, value); + return true; + } + + /* ------------------------------------------------------------------------ */ + virtual void printself(std::ostream & stream) const { + stream << "RemoveDamagedWithDamageRateWeightFunction [damage_limit: " << damage_limit_with_damage_rate << "]"; + } + +private: + /// limit at which a point is considered as complitely broken + Real damage_limit_with_damage_rate; + + /// internal pointer to the current damage vector + const Vector * selected_damage_with_damage_rate; + + /// internal pointer to the current damage rate vector + const Vector * selected_damage_rate_with_damage_rate; +}; /* -------------------------------------------------------------------------- */ /* Stress Based Weight */ /* -------------------------------------------------------------------------- */ template class StressBasedWeightFunction : public BaseWeightFunction { public: StressBasedWeightFunction(Material & material); void init(); virtual void updateInternals(__attribute__((unused)) const ByElementTypeReal & quadrature_points_coordinates) { updatePrincipalStress(_not_ghost); updatePrincipalStress(_ghost); }; void updatePrincipalStress(GhostType ghost_type); inline void updateQuadraturePointsCoordinates(ByElementTypeReal & quadrature_points_coordinates); inline void selectType(ElementType type1, GhostType ghost_type1, ElementType type2, GhostType ghost_type2); inline Real operator()(Real r, QuadraturePoint & q1, QuadraturePoint & q2); inline Real computeRhoSquare(Real r, types::RVector & eigs, types::Matrix & eigenvects, types::RVector & x_s); /* ------------------------------------------------------------------------ */ bool parseParam(const std::string & key, const std::string & value) { std::stringstream sstr(value); if(key == "ft") { sstr >> ft; } else return BaseWeightFunction::parseParam(key, value); return true; } /* ------------------------------------------------------------------------ */ virtual void printself(std::ostream & stream) const { stream << "StressBasedWeightFunction [ft: " << ft << "]"; } private: Real ft; ByElementTypeReal stress_diag; Vector * selected_stress_diag; ByElementTypeReal stress_base; Vector * selected_stress_base; ByElementTypeReal quadrature_points_coordinates; Vector * selected_position_1; Vector * selected_position_2; ByElementTypeReal characteristic_size; Vector * selected_characteristic_size; }; template inline std::ostream & operator <<(std::ostream & stream, const BaseWeightFunction & _this) { _this.printself(stream); return stream; } template inline std::ostream & operator <<(std::ostream & stream, const BaseWeightFunction * _this) { _this->printself(stream); return stream; } template inline std::ostream & operator >>(std::ostream & stream, __attribute__((unused)) const BaseWeightFunction * _this) { AKANTU_DEBUG_TO_IMPLEMENT(); return stream; } template inline std::ostream & operator >>(std::ostream & stream, __attribute__((unused)) const RemoveDamagedWeightFunction * _this) { AKANTU_DEBUG_TO_IMPLEMENT(); return stream; } template inline std::ostream & operator >>(std::ostream & stream, __attribute__((unused)) const DamagedWeightFunction * _this) { AKANTU_DEBUG_TO_IMPLEMENT(); return stream; } template inline std::ostream & operator >>(std::ostream & stream, __attribute__((unused)) const StressBasedWeightFunction * _this) { AKANTU_DEBUG_TO_IMPLEMENT(); return stream; } #include "weight_function_tmpl.hh" __END_AKANTU__ #endif /* __AKANTU_WEIGHT_FUNCTION_HH__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh index 2ef73ff3e..d34c12dff 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,557 +1,557 @@ /** * @file solid_mechanics_model.hh * @author Nicolas Richart * @date[creation] Thu Jul 22 11:51:06 2010 * @date[last modification] Thu Oct 14 14:00:06 2010 * * @brief Model of Solid Mechanics * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_MECHANICS_MODEL_HH__ #define __AKANTU_SOLID_MECHANICS_MODEL_HH__ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "model.hh" #include "data_accessor.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #include "integrator_cohesive.hh" #include "shape_cohesive.hh" #include "aka_types.hh" #include "integration_scheme_2nd_order.hh" #include "solver.hh" /* -------------------------------------------------------------------------- */ namespace akantu { class Material; class IntegrationScheme2ndOrder; class Contact; class SparseMatrix; } __BEGIN_AKANTU__ class SolidMechanicsModel : public Model, public DataAccessor, public MeshEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef FEMTemplate MyFEMType; typedef FEMTemplate< IntegratorCohesive, ShapeCohesive > MyFEMCohesiveType; SolidMechanicsModel(Mesh & mesh, UInt spatial_dimension = 0, const ID & id = "solid_mechanics_model", const MemoryID & memory_id = 0); virtual ~SolidMechanicsModel(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize completely the model virtual void initFull(std::string material_file, AnalysisMethod method = _explicit_dynamic); /// initialize the fem object needed for boundary conditions void initFEMBoundary(bool create_surface = true); /// register the tags associated with the parallel synchronizer void initParallel(MeshPartition * partition, DataAccessor * data_accessor=NULL); /// allocate all vectors void initVectors(); /// initialize all internal arrays for materials void initMaterials(); /// initialize the model virtual void initModel(); /// init PBC synchronizer void initPBC(); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* PBC */ /* ------------------------------------------------------------------------ */ public: /// change the equation number for proper assembly when using PBC void changeEquationNumberforPBC(std::map & pbc_pair); /// synchronize Residual for output void synchronizeResidual(); protected: /// register PBC synchronizer void registerPBCSynchronizer(); /* ------------------------------------------------------------------------ */ /* Explicit */ /* ------------------------------------------------------------------------ */ public: /// initialize the stuff for the explicit scheme void initExplicit(); /// initialize the array needed by updateResidual (residual, current_position) void initializeUpdateResidualData(); /// update the current position vector void updateCurrentPosition(); /// assemble the residual for the explicit scheme void updateResidual(bool need_initialize = true); /** * \brief compute the acceleration from the residual * this function is the explicit equivalent to solveDynamic in implicit * In the case of lumped mass just divide the residual by the mass * In the case of not lumped mass call solveDynamic<_acceleration_corrector> */ void updateAcceleration(); /// Solve the system @f[ A x = \alpha b @f] with A a lumped matrix void solveLumped(Vector & x, const Vector & A, const Vector & b, const Vector & boundary, Real alpha); /// explicit integration predictor void explicitPred(); /// explicit integration corrector void explicitCorr(); /* ------------------------------------------------------------------------ */ /* Implicit */ /* ------------------------------------------------------------------------ */ public: /// initialize the solver and the jacobian_matrix (called by initImplicit) void initSolver(SolverOptions & options = _solver_no_options); /// initialize the stuff for the implicit solver void initImplicit(bool dynamic = false, SolverOptions & solver_options = _solver_no_options); /// solve Ma = f to get the initial acceleration void initialAcceleration(); /// assemble the stiffness matrix void assembleStiffnessMatrix(); /// solve @f[ A\delta u = f_ext - f_int @f] in displacement void solveDynamic(); /// solve Ku = f void solveStatic(); /// test the convergence (norm of increment) bool testConvergenceIncrement(Real tolerance); bool testConvergenceIncrement(Real tolerance, Real & error); /// test the convergence (norm of residual) bool testConvergenceResidual(Real tolerance); bool testConvergenceResidual(Real tolerance, Real & error); /// implicit time integration predictor void implicitPred(); /// implicit time integration corrector void implicitCorr(); protected: /// finish the computation of residual to solve in increment void updateResidualInternal(); /// compute A and solve @f[ A\delta u = f_ext - f_int @f] template void solveDynamic(Vector & increment); /* ------------------------------------------------------------------------ */ /* Boundaries (solid_mechanics_model_boundary.cc) */ /* ------------------------------------------------------------------------ */ public: class SurfaceLoadFunctor { public: virtual void traction(__attribute__ ((unused)) const types::Vector & position, __attribute__ ((unused)) types::Vector & traction, __attribute__ ((unused)) const types::Vector & normal, __attribute__ ((unused)) Surface surface_id) { AKANTU_DEBUG_TO_IMPLEMENT(); } virtual void stress(__attribute__ ((unused)) const types::Vector & position, __attribute__ ((unused)) types::Matrix & stress, __attribute__ ((unused)) const types::Vector & normal, __attribute__ ((unused)) Surface surface_id) { AKANTU_DEBUG_TO_IMPLEMENT(); } }; /// compute force vector from a function(x,y,z) that describe stresses void computeForcesFromFunction(BoundaryFunction in_function, BoundaryFunctionType function_type) __attribute__((deprecated)); template void computeForcesFromFunction(Functor & functor, BoundaryFunctionType function_type); /// integrate a force on the boundary by providing a stress tensor void computeForcesByStressTensor(const Vector & stresses, const ElementType & type, const GhostType & ghost_type); /// integrate a force on the boundary by providing a traction vector void computeForcesByTractionVector(const Vector & tractions, const ElementType & type, const GhostType & ghost_type); /// synchronize the ghost element boundaries values void synchronizeBoundaries(); /* ------------------------------------------------------------------------ */ /* Materials (solid_mechanics_model_material.cc) */ /* ------------------------------------------------------------------------ */ public: /// register a material in the dynamic database Material & registerNewMaterial(const ID & mat_type, const std::string & opt_param = ""); template Material & registerNewCustomMaterial(const ID & mat_type, const std::string & opt_param = ""); /// read the material files to instantiate all the materials void readMaterials(const std::string & filename); /// read a custom material with a keyword and class as template template UInt readCustomMaterial(const std::string & filename, const std::string & keyword); /// Use a UIntData in the mesh to specify the material to use per element void setMaterialIDsFromIntData(const std::string & data_name); protected: // /// read properties part of a material file and create the material // template // Material * readMaterialProperties(std::ifstream & infile, // ID mat_id, // UInt ¤t_line); /* ------------------------------------------------------------------------ */ /* Mass (solid_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the mass matrix void assembleMass(); protected: /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumped(GhostType ghost_type); void assembleMass(GhostType ghost_type); /// fill a vector of rho void computeRho(Vector & rho, ElementType type, GhostType ghost_type); /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: inline virtual UInt getNbDataToPack(const Element & element, SynchronizationTag tag) const; inline virtual UInt getNbDataToUnpack(const Element & element, SynchronizationTag tag) const; inline virtual void packData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag) const; inline virtual void unpackData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag); inline virtual UInt getNbDataToPack(SynchronizationTag tag) const; inline virtual UInt getNbDataToUnpack(SynchronizationTag tag) const; inline virtual void packData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag) const; inline virtual void unpackData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag); /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: virtual void onNodesAdded (const Vector & nodes_list); virtual void onElementsAdded (const Vector & nodes_list); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the current value of the time step AKANTU_GET_MACRO(TimeStep, time_step, Real); /// set the value of the time step AKANTU_SET_MACRO(TimeStep, time_step, Real); /// get the value of the conversion from forces/ mass to acceleration AKANTU_GET_MACRO(F_M2A, f_m2a, Real); /// set the value of the conversion from forces/ mass to acceleration AKANTU_SET_MACRO(F_M2A, f_m2a, Real); /// get the SolidMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement, Vector &); /// get the SolidMechanicsModel::current_position vector \warn only consistent /// after a call to SolidMechanicsModel::updateCurrentPosition AKANTU_GET_MACRO(CurrentPosition, *current_position, const Vector &); /// get the SolidMechanicsModel::increment vector \warn only consistent if /// SolidMechanicsModel::setIncrementFlagOn has been called before AKANTU_GET_MACRO(Increment, *increment, Vector &); /// get the lumped SolidMechanicsModel::mass vector AKANTU_GET_MACRO(Mass, *mass, Vector &); /// get the SolidMechanicsModel::velocity vector AKANTU_GET_MACRO(Velocity, *velocity, Vector &); /// get the SolidMechanicsModel::acceleration vector, updated by /// SolidMechanicsModel::updateAcceleration AKANTU_GET_MACRO(Acceleration, *acceleration, Vector &); /// get the SolidMechanicsModel::force vector (boundary forces) AKANTU_GET_MACRO(Force, *force, Vector &); /// get the SolidMechanicsModel::residual vector, computed by /// SolidMechanicsModel::updateResidual AKANTU_GET_MACRO(Residual, *residual, Vector &); /// get the SolidMechanicsModel::boundary vector AKANTU_GET_MACRO(Boundary, *boundary, Vector &); /// get the SolidMechnicsModel::incrementAcceleration vector AKANTU_GET_MACRO(IncrementAcceleration, *increment_acceleration, Vector &); /// get the value of the SolidMechanicsModel::increment_flag AKANTU_GET_MACRO(IncrementFlag, increment_flag, bool); /// get the SolidMechanicsModel::element_material vector corresponding to a /// given akantu::ElementType AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementMaterial, element_material, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementMaterial, element_material, UInt); /// vectors containing local material element index for each global element index AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementIndexByMaterial, element_index_by_material, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementIndexByMaterial, element_index_by_material, UInt); /// get a particular material inline Material & getMaterial(UInt mat_index); inline const Material & getMaterial(UInt mat_index) const; /// give the number of materials inline UInt getNbMaterials() const { return materials.size(); }; /// give the material internal index from its id Int getInternalIndexFromID(const ID & id) const; /// compute the stable time step Real getStableTimeStep(); /// compute the potential energy Real getPotentialEnergy(); /// compute the kinetic energy Real getKineticEnergy(); /// compute the external work (for impose displacement, the velocity should be given too) Real getExternalWork(); /// get the energies Real getEnergy(std::string id); /// set the Contact object AKANTU_SET_MACRO(Contact, contact, Contact *); /** * @brief set the SolidMechanicsModel::increment_flag to on, the activate the * update of the SolidMechanicsModel::increment vector */ void setIncrementFlagOn(); /// get the stiffness matrix AKANTU_GET_MACRO(StiffnessMatrix, *stiffness_matrix, SparseMatrix &); /// get the mass matrix AKANTU_GET_MACRO(MassMatrix, *mass_matrix, SparseMatrix &); inline FEM & getFEMBoundary(std::string name = ""); /// get integrator AKANTU_GET_MACRO(Integrator, *integrator, const IntegrationScheme2ndOrder &); /// get access to the internal solver AKANTU_GET_MACRO(Solver, *solver, Solver &); protected: /// compute the stable time step Real getStableTimeStep(const GhostType & ghost_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// time step Real time_step; /// conversion coefficient form force/mass to acceleration Real f_m2a; /// displacements array Vector * displacement; /// lumped mass array Vector * mass; /// velocities array Vector * velocity; /// accelerations array Vector * acceleration; /// accelerations array Vector * increment_acceleration; /// forces array Vector * force; /// residuals array Vector * residual; /// boundaries array Vector * boundary; /// array of current position used during update residual Vector * current_position; /// mass matrix SparseMatrix * mass_matrix; /// velocity damping matrix SparseMatrix * velocity_damping_matrix; /// stiffness matrix SparseMatrix * stiffness_matrix; /// jacobian matrix @f[A = cM + dD + K@f] with @f[c = \frac{1}{\beta \Delta t^2}, d = \frac{\gamma}{\beta \Delta t} @f] SparseMatrix * jacobian_matrix; /// materials of all elements ByElementTypeUInt element_material; /// vectors containing local material element index for each global element index ByElementTypeUInt element_index_by_material; /// list of used materials std::vector materials; /// integration scheme of second order used IntegrationScheme2ndOrder * integrator; /// increment of displacement Vector * increment; /// flag defining if the increment must be computed or not bool increment_flag; /// solver for implicit Solver * solver; /// object to resolve the contact Contact * contact; /// the spatial dimension UInt spatial_dimension; /// Mesh Mesh & mesh; AnalysisMethod method; }; __END_AKANTU__ /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ -#include "material.hh" #include "parser.hh" +#include "material.hh" __BEGIN_AKANTU__ #include "solid_mechanics_model_tmpl.hh" #if defined (AKANTU_INCLUDE_INLINE_IMPL) # include "solid_mechanics_model_inline_impl.cc" #endif /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const SolidMechanicsModel & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_SOLID_MECHANICS_MODEL_HH__ */ diff --git a/src/synchronizer/synchronizer.hh b/src/synchronizer/synchronizer.hh index 1afe62b5c..ba96706e9 100644 --- a/src/synchronizer/synchronizer.hh +++ b/src/synchronizer/synchronizer.hh @@ -1,100 +1,101 @@ /** * @file synchronizer.hh * @author Nicolas Richart * @date Mon Aug 23 13:48:37 2010 * * @brief interface for communicator and pbc synchronizers * * @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_SYNCHRONIZER_HH__ #define __AKANTU_SYNCHRONIZER_HH__ /* -------------------------------------------------------------------------- */ #include "aka_memory.hh" #include "data_accessor.hh" /* -------------------------------------------------------------------------- */ namespace akantu { class GhostSynchronizer; } __BEGIN_AKANTU__ class Synchronizer : protected Memory { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Synchronizer(SynchronizerID id = "synchronizer", MemoryID memory_id = 0); virtual ~Synchronizer() { }; - virtual void printself(std::ostream & stream, int indent = 0) const {}; + virtual void printself(__attribute__((unused)) std::ostream & stream, + __attribute__((unused)) int indent = 0) const {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// synchronize ghosts void synchronize(DataAccessor & data_accessor,SynchronizationTag tag); /// asynchronous synchronization of ghosts virtual void asynchronousSynchronize(DataAccessor & data_accessor,SynchronizationTag tag) = 0; /// wait end of asynchronous synchronization of ghosts virtual void waitEndSynchronize(DataAccessor & data_accessor,SynchronizationTag tag) = 0; /// compute buffer size for a given tag and data accessor virtual void computeBufferSize(DataAccessor & data_accessor, SynchronizationTag tag)=0; protected: /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id of the synchronizer SynchronizerID id; }; /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Synchronizer & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_SYNCHRONIZER_HH__ */