diff --git a/src/mesh/mesh.hh b/src/mesh/mesh.hh index c5427a463..547a7e293 100644 --- a/src/mesh/mesh.hh +++ b/src/mesh/mesh.hh @@ -1,701 +1,702 @@ /** * @file mesh.hh * * @author Guillaume Anciaux * @author Dana Christen * @author David Simon Kammer * @author Nicolas Richart * @author Marco Vocialta * * @date creation: Fri Jun 18 2010 * @date last modification: Mon Feb 19 2018 * * @brief the class representing the meshes * * * Copyright (©) 2010-2018 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_array.hh" #include "aka_bbox.hh" #include "aka_event_handler_manager.hh" #include "aka_memory.hh" #include "communicator.hh" #include "dumpable.hh" #include "element.hh" #include "element_class.hh" #include "element_type_map.hh" #include "group_manager.hh" #include "mesh_data.hh" #include "mesh_events.hh" /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ namespace akantu { class ElementSynchronizer; class NodeSynchronizer; class PeriodicNodeSynchronizer; class MeshGlobalDataUpdater; } // namespace akantu namespace akantu { namespace { DECLARE_NAMED_ARGUMENT(communicator); DECLARE_NAMED_ARGUMENT(edge_weight_function); DECLARE_NAMED_ARGUMENT(vertex_weight_function); } // namespace /* -------------------------------------------------------------------------- */ /* Mesh */ /* -------------------------------------------------------------------------- */ /** * @class Mesh mesh.hh * * This class contaisn the coordinates of the nodes in the Mesh.nodes * akant::Array, and the connectivity. The connectivity are stored in by element * types. * * In order to loop on all element you have to loop on all types like this : * @code{.cpp} for(auto & type : mesh.elementTypes()) { UInt nb_element = mesh.getNbElement(type); const Array & conn = mesh.getConnectivity(type); for(UInt e = 0; e < nb_element; ++e) { ... } } or for_each_element(mesh, [](Element & element) { std::cout << element << std::endl }); @endcode */ class Mesh : protected Memory, public EventHandlerManager, public GroupManager, public MeshData, public Dumpable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ private: /// default constructor used for chaining, the last parameter is just to /// differentiate constructors Mesh(UInt spatial_dimension, const ID & id, const MemoryID & memory_id, Communicator & communicator); public: /// constructor that create nodes coordinates array Mesh(UInt spatial_dimension, const ID & id = "mesh", const MemoryID & memory_id = 0); /// mesh not distributed and not using the default communicator Mesh(UInt spatial_dimension, Communicator & communicator, 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, const std::shared_ptr> & nodes, const ID & id = "mesh", const MemoryID & memory_id = 0); ~Mesh() override; /// read the mesh from a file void read(const std::string & filename, const MeshIOType & mesh_io_type = _miot_auto); /// write the mesh to a file void write(const std::string & filename, const MeshIOType & mesh_io_type = _miot_auto); protected: void makeReady(); private: /// initialize the connectivity to NULL and other stuff void init(); /// function that computes the bounding box (fills xmin, xmax) void computeBoundingBox(); /* ------------------------------------------------------------------------ */ /* Distributed memory methods and accessors */ /* ------------------------------------------------------------------------ */ public: protected: /// patitionate the mesh among the processors involved in their computation virtual void distributeImpl( Communicator & communicator, const std::function & edge_weight_function, const std::function & vertex_weight_function); public: /// with the arguments to pass to the partitionner template std::enable_if_t::value> distribute(pack &&... _pack) { distributeImpl( OPTIONAL_NAMED_ARG(communicator, Communicator::getWorldCommunicator()), OPTIONAL_NAMED_ARG(edge_weight_function, [](auto &&, auto &&) { return 1; }), OPTIONAL_NAMED_ARG(vertex_weight_function, [](auto &&) { return 1; })); } /// defines is the mesh is distributed or not inline bool isDistributed() const { return this->is_distributed; } /* ------------------------------------------------------------------------ */ /* Periodicity methods and accessors */ /* ------------------------------------------------------------------------ */ public: /// set the periodicity in a given direction void makePeriodic(const SpatialDirection & direction); void makePeriodic(const SpatialDirection & direction, const ID & list_1, const ID & list_2); protected: void makePeriodic(const SpatialDirection & direction, const Array & list_left, const Array & list_right); /// Removes the face that the mesh is periodic void wipePeriodicInfo(); inline void addPeriodicSlave(UInt slave, UInt master); template void synchronizePeriodicSlaveDataWithMaster(Array & data); // update the periodic synchronizer (creates it if it does not exists) void updatePeriodicSynchronizer(); public: /// defines if the mesh is periodic or not inline bool isPeriodic() const { return this->is_periodic; } inline bool isPeriodic(const SpatialDirection & /*direction*/) const { return this->is_periodic; } class PeriodicSlaves; /// get the master node for a given slave nodes, except if node not a slave inline UInt getPeriodicMaster(UInt slave) const; /// get an iterable list of slaves for a given master node inline decltype(auto) getPeriodicSlaves(UInt master) const; /* ------------------------------------------------------------------------ */ /* General Methods */ /* ------------------------------------------------------------------------ */ public: /// function to print the containt of the class void printself(std::ostream & stream, int indent = 0) const override; /// extract coordinates of nodes from an element template inline void extractNodalValuesFromElement(const Array & nodal_values, T * local_coord, const 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); /// add a Array of connectivity for the given ElementType and GhostType . inline void addConnectivityType(ElementType type, GhostType ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ template inline void sendEvent(Event & event) { // if(event.getList().size() != 0) EventHandlerManager::sendEvent(event); } /// prepare the event to remove the elements listed void eraseElements(const Array & elements); /* ------------------------------------------------------------------------ */ template inline void removeNodesFromArray(Array & vect, const Array & new_numbering); /// initialize normals void initNormals(); /// init facets' mesh Mesh & initMeshFacets(const ID & id = "mesh_facets"); /// define parent mesh void defineMeshParent(const Mesh & mesh); /// get global connectivity array void getGlobalConnectivity(ElementTypeMapArray & global_connectivity); public: void getAssociatedElements(const Array & node_list, Array & elements); private: /// fills the nodes_to_elements structure void fillNodesToElements(); /// update the global ids, nodes type, ... std::tuple updateGlobalData(NewNodesEvent & nodes_event, NewElementsEvent & elements_event); void registerGlobalDataUpdater( std::unique_ptr && global_data_updater); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the id of the mesh AKANTU_GET_MACRO(ID, Memory::id, const ID &); /// get the id of the mesh AKANTU_GET_MACRO(MemoryID, Memory::memory_id, const MemoryID &); /// get the spatial dimension of the mesh = number of component of the /// coordinates AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the nodes Array aka coordinates AKANTU_GET_MACRO(Nodes, *nodes, const Array &); AKANTU_GET_MACRO_NOT_CONST(Nodes, *nodes, Array &); /// get the normals for the elements AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Normals, normals, Real); /// get the number of nodes AKANTU_GET_MACRO(NbNodes, nodes->size(), UInt); /// get the Array of global ids of the nodes (only used in parallel) AKANTU_GET_MACRO(GlobalNodesIds, *nodes_global_ids, const Array &); // AKANTU_GET_MACRO_NOT_CONST(GlobalNodesIds, *nodes_global_ids, Array // &); /// get the global id of a node inline UInt getNodeGlobalId(UInt local_id) const; /// get the global id of a node inline UInt getNodeLocalId(UInt global_id) const; /// get the global number of nodes inline UInt getNbGlobalNodes() const; /// get the nodes type Array AKANTU_GET_MACRO(NodesFlags, *nodes_flags, const Array &); protected: AKANTU_GET_MACRO_NOT_CONST(NodesFlags, *nodes_flags, Array &); public: inline NodeFlag getNodeFlag(UInt local_id) const; inline Int getNodePrank(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; inline bool isPeriodicSlave(UInt n) const; inline bool isPeriodicMaster(UInt n) const; const Vector & getLowerBounds() const { return bbox.getLowerBounds(); } const Vector & getUpperBounds() const { return bbox.getUpperBounds(); } AKANTU_GET_MACRO(BBox, bbox, const BBox &); const Vector & getLocalLowerBounds() const { return bbox_local.getLowerBounds(); } const Vector & getLocalUpperBounds() const { return bbox_local.getUpperBounds(); } AKANTU_GET_MACRO(LocalBBox, bbox_local, const BBox &); /// get the connectivity Array for a given type AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Connectivity, connectivities, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Connectivity, connectivities, UInt); AKANTU_GET_MACRO(Connectivities, connectivities, const ElementTypeMapArray &); /// get the number of element of a type in the mesh inline UInt getNbElement(ElementType type, GhostType ghost_type = _not_ghost) const; /// get the number of element for a given ghost_type and a given dimension inline UInt getNbElement(UInt spatial_dimension = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const; /// compute the barycenter of a given element inline void getBarycenter(const Element & element, Vector & barycenter) const; void getBarycenters(Array & barycenter, ElementType type, GhostType ghost_type) const; /// compute the center of an inscribed circle of a triangular element inline void getIncenter(const Element & element, Vector & incenter) const; /// get the element connected to a subelement (element of lower dimension) const auto & getElementToSubelement() const; /// get the element connected to a subelement const auto & getElementToSubelement(ElementType el_type, GhostType ghost_type = _not_ghost) const; /// get the elements connected to a subelement const auto & getElementToSubelement(const Element & element) const; /// get the subelement (element of lower dimension) connected to a element const auto & getSubelementToElement() const; /// get the subelement connected to an element const auto & getSubelementToElement(ElementType el_type, GhostType ghost_type = _not_ghost) const; /// get the subelement (element of lower dimension) connected to a element VectorProxy getSubelementToElement(const Element & element) const; /// get connectivity of a given element inline VectorProxy getConnectivity(const Element & element) const; inline Vector getConnectivityWithPeriodicity(const Element & element) const; protected: /// get the element connected to a subelement (element of lower dimension) auto & getElementToSubelementNC(); auto & getSubelementToElementNC(); inline auto & getElementToSubelementNC(const Element & element); inline VectorProxy getSubelementToElementNC(const Element & element); /// get the element connected to a subelement auto & getElementToSubelementNC(ElementType el_type, GhostType ghost_type = _not_ghost); /// get the subelement connected to an element auto & getSubelementToElementNC(ElementType el_type, GhostType ghost_type = _not_ghost); inline VectorProxy getConnectivityNC(const Element & element); public: /// get a name field associated to the mesh template inline const Array & getData(const ID & data_name, ElementType el_type, GhostType ghost_type = _not_ghost) const; /// get a name field associated to the mesh template inline Array & getData(const ID & data_name, ElementType el_type, GhostType ghost_type = _not_ghost); /// get a name field associated to the mesh template inline const ElementTypeMapArray & getData(const ID & data_name) const; /// get a name field associated to the mesh template inline ElementTypeMapArray & getData(const ID & data_name); template ElementTypeMap getNbDataPerElem(ElementTypeMapArray & array); template std::shared_ptr createFieldFromAttachedData(const std::string & field_id, const std::string & group_name, ElementKind element_kind); /// templated getter returning the pointer to data in MeshData (modifiable) template inline Array & getDataPointer(const std::string & data_name, ElementType el_type, GhostType ghost_type = _not_ghost, UInt nb_component = 1, bool size_to_nb_element = true, bool resize_with_parent = false); template inline Array & getDataPointer(const ID & data_name, ElementType el_type, GhostType ghost_type, UInt nb_component, bool size_to_nb_element, bool resize_with_parent, const T & defaul_); /// Facets mesh accessor inline bool hasMeshFacets(); inline const Mesh & getMeshFacets() const; inline Mesh & getMeshFacets(); inline auto hasMeshFacets() const { return mesh_facets != nullptr; } /// Parent mesh accessor inline const Mesh & getMeshParent() const; inline bool isMeshFacets() const { return this->is_mesh_facets; } /// return the dumper from a group and and a dumper name DumperIOHelper & getGroupDumper(const std::string & dumper_name, const std::string & group_name); /* ------------------------------------------------------------------------ */ /* Wrappers on ElementClass functions */ /* ------------------------------------------------------------------------ */ public: /// get the number of nodes per element for a given element type static inline UInt getNbNodesPerElement(ElementType type); /// get the number of nodes per element for a given element type considered as /// a first order element static inline ElementType getP1ElementType(ElementType type); /// get the kind of the element type static inline ElementKind getKind(ElementType type); /// get spatial dimension of a type of element static inline UInt getSpatialDimension(ElementType type); /// get number of facets of a given element type static inline UInt getNbFacetsPerElement(ElementType type); /// get number of facets of a given element type static inline UInt getNbFacetsPerElement(ElementType type, UInt t); /// get local connectivity of a facet for a given facet type static inline auto getFacetLocalConnectivity(ElementType type, UInt t = 0); /// get connectivity of facets for a given element inline auto getFacetConnectivity(const Element & element, UInt t = 0) const; /// get the number of type of the surface element associated to a given /// element type static inline UInt getNbFacetTypes(ElementType type, UInt t = 0); /// get the type of the surface element associated to a given element static inline constexpr auto getFacetType(ElementType type, UInt t = 0); /// get all the type of the surface element associated to a given element static inline constexpr auto getAllFacetTypes(ElementType type); /// get the number of nodes in the given element list static inline UInt getNbNodesPerElementList(const Array & elements); /* ------------------------------------------------------------------------ */ /* Element type Iterator */ /* ------------------------------------------------------------------------ */ using type_iterator [[deprecated]] = ElementTypeMapArray::type_iterator; using ElementTypesIteratorHelper = ElementTypeMapArray::ElementTypesIteratorHelper; template ElementTypesIteratorHelper elementTypes(pack &&... _pack) const; [[deprecated("Use elementTypes instead")]] inline decltype(auto) firstType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const { return connectivities.elementTypes(dim, ghost_type, kind).begin(); } [[deprecated("Use elementTypes instead")]] inline decltype(auto) lastType(UInt dim = _all_dimensions, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const { return connectivities.elementTypes(dim, ghost_type, kind).end(); } AKANTU_GET_MACRO(ElementSynchronizer, *element_synchronizer, const ElementSynchronizer &); AKANTU_GET_MACRO_NOT_CONST(ElementSynchronizer, *element_synchronizer, ElementSynchronizer &); AKANTU_GET_MACRO(NodeSynchronizer, *node_synchronizer, const NodeSynchronizer &); AKANTU_GET_MACRO_NOT_CONST(NodeSynchronizer, *node_synchronizer, NodeSynchronizer &); AKANTU_GET_MACRO(PeriodicNodeSynchronizer, *periodic_node_synchronizer, const PeriodicNodeSynchronizer &); AKANTU_GET_MACRO_NOT_CONST(PeriodicNodeSynchronizer, *periodic_node_synchronizer, PeriodicNodeSynchronizer &); // AKANTU_GET_MACRO_NOT_CONST(Communicator, *communicator, StaticCommunicator // &); AKANTU_GET_MACRO(Communicator, *communicator, const auto &); AKANTU_GET_MACRO_NOT_CONST(Communicator, *communicator, auto &); AKANTU_GET_MACRO(PeriodicMasterSlaves, periodic_master_slave, const auto &); + AKANTU_GET_MACRO(PeriodicSlaveMaster, periodic_slave_master, const auto &); /* ------------------------------------------------------------------------ */ /* Private methods for friends */ /* ------------------------------------------------------------------------ */ private: friend class MeshAccessor; friend class MeshUtils; AKANTU_GET_MACRO(NodesPointer, *nodes, Array &); /// get a pointer to the nodes_global_ids Array and create it if /// necessary inline Array & getNodesGlobalIdsPointer(); /// get a pointer to the nodes_type Array and create it if necessary inline Array & getNodesFlagsPointer(); /// get a pointer to the connectivity Array for the given type and create it /// if necessary inline Array & getConnectivityPointer(ElementType type, GhostType ghost_type = _not_ghost); /// get the ghost element counter inline Array & getGhostsCounters(ElementType type, GhostType ghost_type = _ghost) { AKANTU_DEBUG_ASSERT(ghost_type != _not_ghost, "No ghost counter for _not_ghost elements"); return ghosts_counters(type, ghost_type); } /// get a pointer to the element_to_subelement Array for the given type and /// create it if necessary inline Array> & getElementToSubelementPointer(ElementType type, GhostType ghost_type = _not_ghost); /// get a pointer to the subelement_to_element Array for the given type and /// create it if necessary inline Array & getSubelementToElementPointer(ElementType type, GhostType ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// array of the nodes coordinates std::shared_ptr> nodes; /// global node ids std::shared_ptr> nodes_global_ids; /// node flags (shared/periodic/...) std::shared_ptr> nodes_flags; /// processor handling the node when not local or master std::unordered_map nodes_prank; /// global number of nodes; UInt nb_global_nodes{0}; /// all class of elements present in this mesh (for heterogenous meshes) ElementTypeMapArray connectivities; /// count the references on ghost elements ElementTypeMapArray ghosts_counters; /// map to normals for all class of elements present in this mesh ElementTypeMapArray normals; /// the spatial dimension of this mesh UInt spatial_dimension{0}; /// size covered by the mesh on each direction Vector size; /// global bounding box BBox bbox; /// local bounding box BBox bbox_local; /// Extra data loaded from the mesh file // MeshData mesh_data; /// facets' mesh std::unique_ptr mesh_facets; /// parent mesh (this is set for mesh_facets meshes) const Mesh * mesh_parent{nullptr}; /// defines if current mesh is mesh_facets or not bool is_mesh_facets{false}; /// defines if the mesh is centralized or distributed bool is_distributed{false}; /// defines if the mesh is periodic bool is_periodic{false}; /// Communicator on which mesh is distributed Communicator * communicator; /// Element synchronizer std::unique_ptr element_synchronizer; /// Node synchronizer std::unique_ptr node_synchronizer; /// Node synchronizer for periodic nodes std::unique_ptr periodic_node_synchronizer; using NodesToElements = std::vector>>; /// class to update global data using external knowledge std::unique_ptr global_data_updater; /// This info is stored to simplify the dynamic changes NodesToElements nodes_to_elements; /// periodicity local info std::unordered_map periodic_slave_master; std::unordered_multimap periodic_master_slave; }; /// standard output stream operator inline std::ostream & operator<<(std::ostream & stream, const Mesh & _this) { _this.printself(stream); return stream; } } // namespace akantu /* -------------------------------------------------------------------------- */ /* Inline functions */ /* -------------------------------------------------------------------------- */ #include "element_type_map_tmpl.hh" #include "mesh_inline_impl.hh" #endif /* AKANTU_MESH_HH_ */ diff --git a/src/mesh/mesh_accessor.hh b/src/mesh/mesh_accessor.hh index 8e8ecd599..36ba24c88 100644 --- a/src/mesh/mesh_accessor.hh +++ b/src/mesh/mesh_accessor.hh @@ -1,210 +1,214 @@ /** * @file mesh_accessor.hh * * @author Nicolas Richart * * @date creation: Tue Jun 30 2015 * @date last modification: Tue Sep 19 2017 * * @brief this class allow to access some private member of mesh it is used for * IO for examples * * * Copyright (©) 2015-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "mesh.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_MESH_ACCESSOR_HH_ #define AKANTU_MESH_ACCESSOR_HH_ namespace akantu { class NodeSynchronizer; class ElementSynchronizer; class MeshGlobalDataUpdater; } // namespace akantu namespace akantu { class MeshAccessor { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: explicit MeshAccessor(Mesh & mesh) : _mesh(mesh) {} virtual ~MeshAccessor() = default; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the global number of nodes inline UInt getNbGlobalNodes() const { return this->_mesh.nb_global_nodes; } /// set the global number of nodes inline void setNbGlobalNodes(UInt nb_global_nodes) { this->_mesh.nb_global_nodes = nb_global_nodes; } /// set the mesh as being distributed inline void setDistributed() { this->_mesh.is_distributed = true; } /// get a pointer to the nodes_global_ids Array and create it if /// necessary inline auto & getNodesGlobalIds() { return this->_mesh.getNodesGlobalIdsPointer(); } /// get a pointer to the nodes_type Array and create it if necessary inline auto & getNodesFlags() { return this->_mesh.getNodesFlags(); } /// get a pointer to the nodes_type Array and create it if necessary inline void setNodePrank(UInt node, Int prank) { this->_mesh.nodes_prank[node] = prank; } /// get a pointer to the coordinates Array inline auto & getNodes() { return this->_mesh.getNodesPointer(); } /// get a pointer to the coordinates Array inline auto getNodesSharedPtr() { return this->_mesh.nodes; } /// get the connectivities inline auto & getConnectivities() { return this->_mesh.connectivities; } /// get the connectivity Array for the given type and create it /// if necessary inline auto & getConnectivity(ElementType type, GhostType ghost_type = _not_ghost) { return this->_mesh.getConnectivityPointer(type, ghost_type); } /// get the connectivity for the given element inline decltype(auto) getConnectivity(const Element & element) { return this->_mesh.getConnectivityNC(element); } /// get the ghost element counter inline auto & getGhostsCounters(ElementType type, GhostType ghost_type = _ghost) { return this->_mesh.getGhostsCounters(type, ghost_type); } /// get the element_to_subelement Array for the given type and /// create it if necessary inline auto & getElementToSubelement(ElementType type, GhostType ghost_type = _not_ghost) { return this->_mesh.getElementToSubelementPointer(type, ghost_type); } inline decltype(auto) getElementToSubelementNC(const ElementType & type, const GhostType & ghost_type = _not_ghost) { return this->_mesh.getElementToSubelementNC(type, ghost_type); } /// get the subelement_to_element Array for the given type and /// create it if necessary inline auto & getSubelementToElement(ElementType type, GhostType ghost_type = _not_ghost) { return this->_mesh.getSubelementToElementPointer(type, ghost_type); } inline decltype(auto) getSubelementToElementNC(const ElementType & type, const GhostType & ghost_type = _not_ghost) { return this->_mesh.getSubelementToElementNC(type, ghost_type); } /// get the element_to_subelement, creates it if necessary inline decltype(auto) getElementToSubelement() { return this->_mesh.getElementToSubelementNC(); } /// get subelement_to_element, creates it if necessary inline decltype(auto) getSubelementToElement() { return this->_mesh.getSubelementToElementNC(); } /// get a pointer to the element_to_subelement Array for element and /// create it if necessary inline decltype(auto) getElementToSubelement(const Element & element) { return this->_mesh.getElementToSubelementNC(element); } /// get a pointer to the subelement_to_element Array for the given element and /// create it if necessary inline decltype(auto) getSubelementToElement(const Element & element) { return this->_mesh.getSubelementToElementNC(element); } template inline auto & getData(const std::string & data_name, ElementType el_type, GhostType ghost_type = _not_ghost, UInt nb_component = 1, bool size_to_nb_element = true, bool resize_with_parent = false) { return this->_mesh.getDataPointer(data_name, el_type, ghost_type, nb_component, size_to_nb_element, resize_with_parent); } /// get the node synchonizer auto & getNodeSynchronizer() { return *this->_mesh.node_synchronizer; } /// get the element synchonizer auto & getElementSynchronizer() { return *this->_mesh.element_synchronizer; } decltype(auto) updateGlobalData(NewNodesEvent & nodes_event, NewElementsEvent & elements_event) { return this->_mesh.updateGlobalData(nodes_event, elements_event); } void registerGlobalDataUpdater( std::unique_ptr && global_data_updater) { this->_mesh.registerGlobalDataUpdater( std::forward>( global_data_updater)); } /* ------------------------------------------------------------------------ */ void makeReady() { this->_mesh.makeReady(); } /* ------------------------------------------------------------------------ */ void addPeriodicSlave(UInt slave, UInt master) { this->_mesh.addPeriodicSlave(slave, master); } void markMeshPeriodic() { for (UInt s : arange(this->_mesh.spatial_dimension)) { this->_mesh.is_periodic |= 1 << s; } } + void updatePeriodicSynchronizer() { + this->_mesh.updatePeriodicSynchronizer(); + } + void wipePeriodicInfo() { this->_mesh.wipePeriodicInfo(); } private: Mesh & _mesh; }; } // namespace akantu #endif /* AKANTU_MESH_ACCESSOR_HH_ */ diff --git a/src/model/common/dof_manager/dof_manager.cc b/src/model/common/dof_manager/dof_manager.cc index 0fc775deb..66d90f29a 100644 --- a/src/model/common/dof_manager/dof_manager.cc +++ b/src/model/common/dof_manager/dof_manager.cc @@ -1,1018 +1,1030 @@ /** * @file dof_manager.cc * * @author Nicolas Richart * * @date creation: Tue Aug 18 2015 * @date last modification: Wed Feb 21 2018 * * @brief Implementation of the common parts of the DOFManagers * * * Copyright (©) 2015-2018 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 "dof_manager.hh" #include "communicator.hh" #include "mesh.hh" #include "mesh_utils.hh" #include "node_group.hh" #include "node_synchronizer.hh" #include "non_linear_solver.hh" #include "periodic_node_synchronizer.hh" #include "time_step_solver.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ DOFManager::DOFManager(const ID & id, const MemoryID & memory_id) : Memory(id, memory_id), dofs_flag(0, 1, std::string(id + ":dofs_type")), global_equation_number(0, 1, "global_equation_number"), communicator(Communicator::getSelfCommunicator()) {} /* -------------------------------------------------------------------------- */ DOFManager::DOFManager(Mesh & mesh, const ID & id, const MemoryID & memory_id) : Memory(id, memory_id), mesh(&mesh), dofs_flag(0, 1, std::string(id + ":dofs_type")), global_equation_number(0, 1, "global_equation_number"), communicator(mesh.getCommunicator()) { this->mesh->registerEventHandler(*this, _ehp_dof_manager); } /* -------------------------------------------------------------------------- */ DOFManager::~DOFManager() = default; /* -------------------------------------------------------------------------- */ std::vector DOFManager::getDOFIDs() const { std::vector keys; for (const auto & dof_data : this->dofs) { keys.push_back(dof_data.first); } return keys; } /* -------------------------------------------------------------------------- */ void DOFManager::assembleElementalArrayLocalArray( const Array & elementary_vect, Array & array_assembeled, ElementType type, GhostType ghost_type, Real scale_factor, const Array & filter_elements) { AKANTU_DEBUG_IN(); UInt nb_element; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_degree_of_freedom = elementary_vect.getNbComponent() / nb_nodes_per_element; UInt * filter_it = nullptr; if (filter_elements != empty_filter) { nb_element = filter_elements.size(); filter_it = filter_elements.storage(); } else { nb_element = this->mesh->getNbElement(type, ghost_type); } AKANTU_DEBUG_ASSERT(elementary_vect.size() == nb_element, "The vector elementary_vect(" << elementary_vect.getID() << ") has not the good size."); const Array & connectivity = this->mesh->getConnectivity(type, ghost_type); Array::const_matrix_iterator elem_it = elementary_vect.begin(nb_degree_of_freedom, nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++elem_it) { UInt element = el; if (filter_it != nullptr) { // conn_it = conn_begin + *filter_it; element = *filter_it; } // const Vector & conn = *conn_it; const Matrix & elemental_val = *elem_it; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt offset_node = connectivity(element, n) * nb_degree_of_freedom; Vector assemble(array_assembeled.storage() + offset_node, nb_degree_of_freedom); Vector elem_val = elemental_val(n); assemble.aXplusY(elem_val, scale_factor); } if (filter_it != nullptr) { ++filter_it; } // else // ++conn_it; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleElementalArrayToResidual( const ID & dof_id, const Array & elementary_vect, ElementType type, GhostType ghost_type, Real scale_factor, const Array & filter_elements) { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_degree_of_freedom = elementary_vect.getNbComponent() / nb_nodes_per_element; Array array_localy_assembeled(this->mesh->getNbNodes(), nb_degree_of_freedom); array_localy_assembeled.zero(); this->assembleElementalArrayLocalArray( elementary_vect, array_localy_assembeled, type, ghost_type, scale_factor, filter_elements); this->assembleToResidual(dof_id, array_localy_assembeled, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleElementalArrayToLumpedMatrix( const ID & dof_id, const Array & elementary_vect, const ID & lumped_mtx, ElementType type, GhostType ghost_type, Real scale_factor, const Array & filter_elements) { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_degree_of_freedom = elementary_vect.getNbComponent() / nb_nodes_per_element; Array array_localy_assembeled(this->mesh->getNbNodes(), nb_degree_of_freedom); array_localy_assembeled.zero(); this->assembleElementalArrayLocalArray( elementary_vect, array_localy_assembeled, type, ghost_type, scale_factor, filter_elements); this->assembleToLumpedMatrix(dof_id, array_localy_assembeled, lumped_mtx, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleMatMulDOFsToResidual(const ID & A_id, Real scale_factor) { for (auto & pair : this->dofs) { const auto & dof_id = pair.first; auto & dof_data = *pair.second; this->assembleMatMulVectToResidual(dof_id, A_id, *dof_data.dof, scale_factor); } } /* -------------------------------------------------------------------------- */ void DOFManager::splitSolutionPerDOFs() { for (auto && data : this->dofs) { auto & dof_data = *data.second; dof_data.solution.resize(dof_data.dof->size() * dof_data.dof->getNbComponent()); this->getSolutionPerDOFs(data.first, dof_data.solution); } } /* -------------------------------------------------------------------------- */ void DOFManager::getSolutionPerDOFs(const ID & dof_id, Array & solution_array) { AKANTU_DEBUG_IN(); this->getArrayPerDOFs(dof_id, this->getSolution(), solution_array); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::getLumpedMatrixPerDOFs(const ID & dof_id, const ID & lumped_mtx, Array & lumped) { AKANTU_DEBUG_IN(); this->getArrayPerDOFs(dof_id, this->getLumpedMatrix(lumped_mtx), lumped); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleToResidual(const ID & dof_id, Array & array_to_assemble, Real scale_factor) { AKANTU_DEBUG_IN(); // this->makeConsistentForPeriodicity(dof_id, array_to_assemble); this->assembleToGlobalArray(dof_id, array_to_assemble, this->getResidual(), scale_factor); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleToLumpedMatrix(const ID & dof_id, Array & array_to_assemble, const ID & lumped_mtx, Real scale_factor) { AKANTU_DEBUG_IN(); // this->makeConsistentForPeriodicity(dof_id, array_to_assemble); auto & lumped = this->getLumpedMatrix(lumped_mtx); this->assembleToGlobalArray(dof_id, array_to_assemble, lumped, scale_factor); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ DOFManager::DOFData::DOFData(const ID & dof_id) : support_type(_dst_generic), group_support("__mesh__"), solution(0, 1, dof_id + ":solution"), local_equation_number(0, 1, dof_id + ":local_equation_number"), associated_nodes(0, 1, dof_id + "associated_nodes") {} /* -------------------------------------------------------------------------- */ DOFManager::DOFData::~DOFData() = default; /* -------------------------------------------------------------------------- */ template auto DOFManager::countDOFsForNodes(const DOFData & dof_data, UInt nb_nodes, Func && getNode) { auto nb_local_dofs = nb_nodes; decltype(nb_local_dofs) nb_pure_local = 0; for (auto n : arange(nb_nodes)) { UInt node = getNode(n); // http://www.open-std.org/jtc1/sc22/open/n2356/conv.html // bool are by convention casted to 0 and 1 when promoted to int nb_pure_local += this->mesh->isLocalOrMasterNode(node); nb_local_dofs -= this->mesh->isPeriodicSlave(node); } const auto & dofs_array = *dof_data.dof; nb_pure_local *= dofs_array.getNbComponent(); nb_local_dofs *= dofs_array.getNbComponent(); return std::make_pair(nb_local_dofs, nb_pure_local); } /* -------------------------------------------------------------------------- */ auto DOFManager::getNewDOFDataInternal(const ID & dof_id) -> DOFData & { auto it = this->dofs.find(dof_id); if (it != this->dofs.end()) { AKANTU_EXCEPTION("This dof array has already been registered"); } std::unique_ptr dof_data_ptr = this->getNewDOFData(dof_id); DOFData & dof_data = *dof_data_ptr; this->dofs[dof_id] = std::move(dof_data_ptr); return dof_data; } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFs(const ID & dof_id, Array & dofs_array, const DOFSupportType & support_type) { auto & dofs_storage = this->getNewDOFDataInternal(dof_id); dofs_storage.support_type = support_type; this->registerDOFsInternal(dof_id, dofs_array); resizeGlobalArrays(); } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFs(const ID & dof_id, Array & dofs_array, const ID & support_group) { auto & dofs_storage = this->getNewDOFDataInternal(dof_id); dofs_storage.support_type = _dst_nodal; dofs_storage.group_support = support_group; this->registerDOFsInternal(dof_id, dofs_array); resizeGlobalArrays(); } /* -------------------------------------------------------------------------- */ std::tuple DOFManager::registerDOFsInternal(const ID & dof_id, Array & dofs_array) { DOFData & dof_data = this->getDOFData(dof_id); dof_data.dof = &dofs_array; UInt nb_local_dofs = 0; UInt nb_pure_local = 0; const auto & support_type = dof_data.support_type; switch (support_type) { case _dst_nodal: { const auto & group = dof_data.group_support; std::function getNode; if (group == "__mesh__") { AKANTU_DEBUG_ASSERT( dofs_array.size() == this->mesh->getNbNodes(), "The array of dof is too short to be associated to nodes."); std::tie(nb_local_dofs, nb_pure_local) = countDOFsForNodes( dof_data, this->mesh->getNbNodes(), [](auto && n) { return n; }); } else { const auto & node_group = this->mesh->getElementGroup(group).getNodeGroup().getNodes(); AKANTU_DEBUG_ASSERT( dofs_array.size() == node_group.size(), "The array of dof is too shot to be associated to nodes."); std::tie(nb_local_dofs, nb_pure_local) = countDOFsForNodes(dof_data, node_group.size(), [&node_group](auto && n) { return node_group(n); }); } break; } case _dst_generic: { nb_local_dofs = nb_pure_local = dofs_array.size() * dofs_array.getNbComponent(); break; } default: { AKANTU_EXCEPTION("This type of dofs is not handled yet."); } } dof_data.local_nb_dofs = nb_local_dofs; dof_data.pure_local_nb_dofs = nb_pure_local; dof_data.ghosts_nb_dofs = nb_local_dofs - nb_pure_local; this->pure_local_system_size += nb_pure_local; this->local_system_size += nb_local_dofs; auto nb_total_pure_local = nb_pure_local; communicator.allReduce(nb_total_pure_local, SynchronizerOperation::_sum); this->system_size += nb_total_pure_local; // updating the dofs data after counting is finished switch (support_type) { case _dst_nodal: { const auto & group = dof_data.group_support; if (group != "__mesh__") { auto & support_nodes = this->mesh->getElementGroup(group).getNodeGroup().getNodes(); this->updateDOFsData( dof_data, nb_local_dofs, nb_pure_local, support_nodes.size(), [&support_nodes](UInt node) -> UInt { return support_nodes[node]; }); } else { this->updateDOFsData(dof_data, nb_local_dofs, nb_pure_local, mesh->getNbNodes(), [](UInt node) -> UInt { return node; }); } break; } case _dst_generic: { this->updateDOFsData(dof_data, nb_local_dofs, nb_pure_local); break; } } return std::make_tuple(nb_local_dofs, nb_pure_local, nb_total_pure_local); } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFsPrevious(const ID & dof_id, Array & array) { DOFData & dof = this->getDOFData(dof_id); if (dof.previous != nullptr) { AKANTU_EXCEPTION("The previous dofs array for " << dof_id << " has already been registered"); } dof.previous = &array; } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFsIncrement(const ID & dof_id, Array & array) { DOFData & dof = this->getDOFData(dof_id); if (dof.increment != nullptr) { AKANTU_EXCEPTION("The dofs increment array for " << dof_id << " has already been registered"); } dof.increment = &array; } /* -------------------------------------------------------------------------- */ void DOFManager::registerDOFsDerivative(const ID & dof_id, UInt order, Array & dofs_derivative) { DOFData & dof = this->getDOFData(dof_id); std::vector *> & derivatives = dof.dof_derivatives; if (derivatives.size() < order) { derivatives.resize(order, nullptr); } else { if (derivatives[order - 1] != nullptr) { AKANTU_EXCEPTION("The dof derivatives of order " << order << " already been registered for this dof (" << dof_id << ")"); } } derivatives[order - 1] = &dofs_derivative; } /* -------------------------------------------------------------------------- */ void DOFManager::registerBlockedDOFs(const ID & dof_id, Array & blocked_dofs) { DOFData & dof = this->getDOFData(dof_id); if (dof.blocked_dofs != nullptr) { AKANTU_EXCEPTION("The blocked dofs array for " << dof_id << " has already been registered"); } dof.blocked_dofs = &blocked_dofs; } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManager::registerSparseMatrix(const ID & matrix_id, std::unique_ptr & matrix) { auto it = this->matrices.find(matrix_id); if (it != this->matrices.end()) { AKANTU_EXCEPTION("The matrix " << matrix_id << " already exists in " << this->id); } auto & ret = *matrix; this->matrices[matrix_id] = std::move(matrix); return ret; } /* -------------------------------------------------------------------------- */ /// Get an instance of a new SparseMatrix SolverVector & DOFManager::registerLumpedMatrix(const ID & matrix_id, std::unique_ptr & matrix) { auto it = this->lumped_matrices.find(matrix_id); if (it != this->lumped_matrices.end()) { AKANTU_EXCEPTION("The lumped matrix " << matrix_id << " already exists in " << this->id); } auto & ret = *matrix; this->lumped_matrices[matrix_id] = std::move(matrix); ret.resize(); return ret; } /* -------------------------------------------------------------------------- */ NonLinearSolver & DOFManager::registerNonLinearSolver( const ID & non_linear_solver_id, std::unique_ptr & non_linear_solver) { NonLinearSolversMap::const_iterator it = this->non_linear_solvers.find(non_linear_solver_id); if (it != this->non_linear_solvers.end()) { AKANTU_EXCEPTION("The non linear solver " << non_linear_solver_id << " already exists in " << this->id); } NonLinearSolver & ret = *non_linear_solver; this->non_linear_solvers[non_linear_solver_id] = std::move(non_linear_solver); return ret; } /* -------------------------------------------------------------------------- */ TimeStepSolver & DOFManager::registerTimeStepSolver( const ID & time_step_solver_id, std::unique_ptr & time_step_solver) { TimeStepSolversMap::const_iterator it = this->time_step_solvers.find(time_step_solver_id); if (it != this->time_step_solvers.end()) { AKANTU_EXCEPTION("The non linear solver " << time_step_solver_id << " already exists in " << this->id); } TimeStepSolver & ret = *time_step_solver; this->time_step_solvers[time_step_solver_id] = std::move(time_step_solver); return ret; } /* -------------------------------------------------------------------------- */ SparseMatrix & DOFManager::getMatrix(const ID & id) { ID matrix_id = this->id + ":mtx:" + id; SparseMatricesMap::const_iterator it = this->matrices.find(matrix_id); if (it == this->matrices.end()) { AKANTU_SILENT_EXCEPTION("The matrix " << matrix_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ bool DOFManager::hasMatrix(const ID & id) const { ID mtx_id = this->id + ":mtx:" + id; auto it = this->matrices.find(mtx_id); return it != this->matrices.end(); } /* -------------------------------------------------------------------------- */ SolverVector & DOFManager::getLumpedMatrix(const ID & id) { ID matrix_id = this->id + ":lumped_mtx:" + id; LumpedMatricesMap::const_iterator it = this->lumped_matrices.find(matrix_id); if (it == this->lumped_matrices.end()) { AKANTU_SILENT_EXCEPTION("The lumped matrix " << matrix_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ const SolverVector & DOFManager::getLumpedMatrix(const ID & id) const { ID matrix_id = this->id + ":lumped_mtx:" + id; auto it = this->lumped_matrices.find(matrix_id); if (it == this->lumped_matrices.end()) { AKANTU_SILENT_EXCEPTION("The lumped matrix " << matrix_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ bool DOFManager::hasLumpedMatrix(const ID & id) const { ID mtx_id = this->id + ":lumped_mtx:" + id; auto it = this->lumped_matrices.find(mtx_id); return it != this->lumped_matrices.end(); } /* -------------------------------------------------------------------------- */ NonLinearSolver & DOFManager::getNonLinearSolver(const ID & id) { ID non_linear_solver_id = this->id + ":nls:" + id; NonLinearSolversMap::const_iterator it = this->non_linear_solvers.find(non_linear_solver_id); if (it == this->non_linear_solvers.end()) { AKANTU_EXCEPTION("The non linear solver " << non_linear_solver_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ bool DOFManager::hasNonLinearSolver(const ID & id) const { ID solver_id = this->id + ":nls:" + id; auto it = this->non_linear_solvers.find(solver_id); return it != this->non_linear_solvers.end(); } /* -------------------------------------------------------------------------- */ TimeStepSolver & DOFManager::getTimeStepSolver(const ID & id) { ID time_step_solver_id = this->id + ":tss:" + id; TimeStepSolversMap::const_iterator it = this->time_step_solvers.find(time_step_solver_id); if (it == this->time_step_solvers.end()) { AKANTU_EXCEPTION("The non linear solver " << time_step_solver_id << " does not exists in " << this->id); } return *(it->second); } /* -------------------------------------------------------------------------- */ bool DOFManager::hasTimeStepSolver(const ID & solver_id) const { ID time_step_solver_id = this->id + ":tss:" + solver_id; auto it = this->time_step_solvers.find(time_step_solver_id); return it != this->time_step_solvers.end(); } /* -------------------------------------------------------------------------- */ void DOFManager::savePreviousDOFs(const ID & dofs_id) { this->getPreviousDOFs(dofs_id).copy(this->getDOFs(dofs_id)); } /* -------------------------------------------------------------------------- */ void DOFManager::zeroResidual() { this->residual->zero(); } /* -------------------------------------------------------------------------- */ void DOFManager::zeroMatrix(const ID & mtx) { this->getMatrix(mtx).zero(); } /* -------------------------------------------------------------------------- */ void DOFManager::zeroLumpedMatrix(const ID & mtx) { this->getLumpedMatrix(mtx).zero(); } /* -------------------------------------------------------------------------- */ /* Mesh Events */ /* -------------------------------------------------------------------------- */ std::pair DOFManager::updateNodalDOFs(const ID & dof_id, const Array & nodes_list) { auto & dof_data = this->getDOFData(dof_id); UInt nb_new_local_dofs; UInt nb_new_pure_local; std::tie(nb_new_local_dofs, nb_new_pure_local) = countDOFsForNodes(dof_data, nodes_list.size(), [&nodes_list](auto && n) { return nodes_list(n); }); this->pure_local_system_size += nb_new_pure_local; this->local_system_size += nb_new_local_dofs; UInt nb_new_global = nb_new_pure_local; communicator.allReduce(nb_new_global, SynchronizerOperation::_sum); this->system_size += nb_new_global; dof_data.solution.resize(local_system_size); updateDOFsData(dof_data, nb_new_local_dofs, nb_new_pure_local, nodes_list.size(), [&nodes_list](UInt pos) -> UInt { return nodes_list[pos]; }); return std::make_pair(nb_new_local_dofs, nb_new_pure_local); } /* -------------------------------------------------------------------------- */ void DOFManager::resizeGlobalArrays() { // resize all relevant arrays this->residual->resize(); this->solution->resize(); this->data_cache->resize(); for (auto & lumped_matrix : lumped_matrices) { lumped_matrix.second->resize(); } for (auto & matrix : matrices) { matrix.second->clearProfile(); matrix.second->resize(this->system_size); } } /* -------------------------------------------------------------------------- */ void DOFManager::onNodesAdded(const Array & nodes_list, const NewNodesEvent & /*unused*/) { for (auto & pair : this->dofs) { const auto & dof_id = pair.first; auto & dof_data = this->getDOFData(dof_id); if (dof_data.support_type != _dst_nodal) { continue; } const auto & group = dof_data.group_support; if (group == "__mesh__") { this->updateNodalDOFs(dof_id, nodes_list); } else { const auto & node_group = this->mesh->getElementGroup(group).getNodeGroup(); Array new_nodes_list; for (const auto & node : nodes_list) { if (node_group.find(node) != UInt(-1)) { new_nodes_list.push_back(node); } } this->updateNodalDOFs(dof_id, new_nodes_list); } } this->resizeGlobalArrays(); } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ class GlobalDOFInfoDataAccessor : public DataAccessor { public: using size_type = typename std::unordered_map>::size_type; GlobalDOFInfoDataAccessor(DOFManager::DOFData & dof_data, DOFManager & dof_manager) : dof_data(dof_data), dof_manager(dof_manager) { for (auto && pair : zip(dof_data.local_equation_number, dof_data.associated_nodes)) { UInt node; Int dof; std::tie(dof, node) = pair; dofs_per_node[node].push_back(dof); } } UInt getNbData(const Array & nodes, const SynchronizationTag & tag) const override { if (tag == SynchronizationTag::_ask_nodes or tag == SynchronizationTag::_giu_global_conn) { return nodes.size() * dof_data.dof->getNbComponent() * sizeof(Int); } return 0; } void packData(CommunicationBuffer & buffer, const Array & nodes, const SynchronizationTag & tag) const override { if (tag == SynchronizationTag::_ask_nodes or tag == SynchronizationTag::_giu_global_conn) { for (const auto & node : nodes) { const auto & dofs = dofs_per_node.at(node); for (const auto & dof : dofs) { buffer << dof_manager.global_equation_number(dof); } } } } void unpackData(CommunicationBuffer & buffer, const Array & nodes, const SynchronizationTag & tag) override { if (tag == SynchronizationTag::_ask_nodes or tag == SynchronizationTag::_giu_global_conn) { for (const auto & node : nodes) { const auto & dofs = dofs_per_node[node]; for (const auto & dof : dofs) { Int global_dof; buffer >> global_dof; AKANTU_DEBUG_ASSERT( (dof_manager.global_equation_number(dof) == -1 or dof_manager.global_equation_number(dof) == global_dof), "This dof already had a global_dof_id which is different from " "the received one. " << dof_manager.global_equation_number(dof) << " != " << global_dof); dof_manager.global_equation_number(dof) = global_dof; dof_manager.global_to_local_mapping[global_dof] = dof; } } } } protected: std::unordered_map> dofs_per_node; DOFManager::DOFData & dof_data; DOFManager & dof_manager; }; /* -------------------------------------------------------------------------- */ auto DOFManager::computeFirstDOFIDs(UInt nb_new_local_dofs, UInt nb_new_pure_local) { // determine the first local/global dof id to use UInt offset = 0; this->communicator.exclusiveScan(nb_new_pure_local, offset); auto first_global_dof_id = this->first_global_dof_id + offset; auto first_local_dof_id = this->local_system_size - nb_new_local_dofs; offset = nb_new_pure_local; this->communicator.allReduce(offset); this->first_global_dof_id += offset; return std::make_pair(first_local_dof_id, first_global_dof_id); } /* -------------------------------------------------------------------------- */ void DOFManager::updateDOFsData(DOFData & dof_data, UInt nb_new_local_dofs, UInt nb_new_pure_local, UInt nb_node, const std::function & getNode) { auto nb_local_dofs_added = nb_node * dof_data.dof->getNbComponent(); auto first_dof_pos = dof_data.local_equation_number.size(); dof_data.local_equation_number.reserve(dof_data.local_equation_number.size() + nb_local_dofs_added); dof_data.associated_nodes.reserve(dof_data.associated_nodes.size() + nb_local_dofs_added); this->dofs_flag.resize(this->local_system_size, NodeFlag::_normal); this->global_equation_number.resize(this->local_system_size, -1); std::unordered_map, UInt> masters_dofs; + // fill in masters_dofs in case masters are already existing nodes + for (auto d : arange(first_dof_pos)) { + auto node = dof_data.associated_nodes(d); + auto is_periodic_master = this->mesh->isPeriodicMaster(node); + if (is_periodic_master) { + auto dof = d % dof_data.dof->getNbComponent(); + auto local_eq_num = dof_data.local_equation_number(d); + masters_dofs.insert( + std::make_pair(std::make_pair(node, dof), local_eq_num)); + } + } + // update per dof info UInt local_eq_num; UInt first_global_dof_id; std::tie(local_eq_num, first_global_dof_id) = computeFirstDOFIDs(nb_new_local_dofs, nb_new_pure_local); for (auto d : arange(nb_local_dofs_added)) { auto node = getNode(d / dof_data.dof->getNbComponent()); auto dof_flag = this->mesh->getNodeFlag(node); dof_data.associated_nodes.push_back(node); auto is_local_dof = this->mesh->isLocalOrMasterNode(node); auto is_periodic_slave = this->mesh->isPeriodicSlave(node); auto is_periodic_master = this->mesh->isPeriodicMaster(node); if (is_periodic_slave) { dof_data.local_equation_number.push_back(-1); continue; } // update equation numbers this->dofs_flag(local_eq_num) = dof_flag; dof_data.local_equation_number.push_back(local_eq_num); if (is_local_dof) { this->global_equation_number(local_eq_num) = first_global_dof_id; this->global_to_local_mapping[first_global_dof_id] = local_eq_num; ++first_global_dof_id; } else { this->global_equation_number(local_eq_num) = -1; } if (is_periodic_master) { auto node = getNode(d / dof_data.dof->getNbComponent()); auto dof = d % dof_data.dof->getNbComponent(); masters_dofs.insert( std::make_pair(std::make_pair(node, dof), local_eq_num)); } ++local_eq_num; } // correct periodic slave equation numbers if (this->mesh->isPeriodic()) { auto assoc_begin = dof_data.associated_nodes.begin(); for (auto d : arange(nb_local_dofs_added)) { auto node = dof_data.associated_nodes(first_dof_pos + d); if (not this->mesh->isPeriodicSlave(node)) { continue; } auto master_node = this->mesh->getPeriodicMaster(node); auto dof = d % dof_data.dof->getNbComponent(); dof_data.local_equation_number(first_dof_pos + d) = - masters_dofs[std::make_pair(master_node, dof)]; + masters_dofs.at(std::make_pair(master_node, dof)); } } // synchronize the global numbering for slaves nodes if (this->mesh->isDistributed()) { GlobalDOFInfoDataAccessor data_accessor(dof_data, *this); if (this->mesh->isPeriodic()) { mesh->getPeriodicNodeSynchronizer().synchronizeOnce( data_accessor, SynchronizationTag::_giu_global_conn); } auto & node_synchronizer = this->mesh->getNodeSynchronizer(); node_synchronizer.synchronizeOnce(data_accessor, SynchronizationTag::_ask_nodes); } } /* -------------------------------------------------------------------------- */ void DOFManager::updateDOFsData(DOFData & dof_data, UInt nb_new_local_dofs, UInt nb_new_pure_local) { dof_data.local_equation_number.reserve(dof_data.local_equation_number.size() + nb_new_local_dofs); UInt first_local_dof_id; UInt first_global_dof_id; std::tie(first_local_dof_id, first_global_dof_id) = computeFirstDOFIDs(nb_new_local_dofs, nb_new_pure_local); this->dofs_flag.resize(this->local_system_size, NodeFlag::_normal); this->global_equation_number.resize(this->local_system_size, -1); // update per dof info for (auto _[[gnu::unused]] : arange(nb_new_local_dofs)) { // update equation numbers this->dofs_flag(first_local_dof_id) = NodeFlag::_normal; dof_data.local_equation_number.push_back(first_local_dof_id); this->global_equation_number(first_local_dof_id) = first_global_dof_id; this->global_to_local_mapping[first_global_dof_id] = first_local_dof_id; ++first_global_dof_id; ++first_local_dof_id; } } /* -------------------------------------------------------------------------- */ void DOFManager::onNodesRemoved(const Array & /*unused*/, const Array & /*unused*/, const RemovedNodesEvent & /*unused*/) {} /* -------------------------------------------------------------------------- */ void DOFManager::onElementsAdded(const Array & /*unused*/, const NewElementsEvent & /*unused*/) {} /* -------------------------------------------------------------------------- */ void DOFManager::onElementsRemoved(const Array & /*unused*/, const ElementTypeMapArray & /*unused*/, const RemovedElementsEvent & /*unused*/) {} /* -------------------------------------------------------------------------- */ void DOFManager::onElementsChanged(const Array & /*unused*/, const Array & /*unused*/, const ElementTypeMapArray & /*unused*/, const ChangedElementsEvent & /*unused*/) {} /* -------------------------------------------------------------------------- */ void DOFManager::updateGlobalBlockedDofs() { this->previous_global_blocked_dofs.copy(this->global_blocked_dofs); this->global_blocked_dofs.reserve(this->local_system_size, 0); this->previous_global_blocked_dofs_release = this->global_blocked_dofs_release; for (auto & pair : dofs) { if (!this->hasBlockedDOFs(pair.first)) { continue; } DOFData & dof_data = *pair.second; for (auto && data : zip(dof_data.getLocalEquationsNumbers(), make_view(*dof_data.blocked_dofs))) { const auto & dof = std::get<0>(data); const auto & is_blocked = std::get<1>(data); if (is_blocked) { this->global_blocked_dofs.push_back(dof); } } } std::sort(this->global_blocked_dofs.begin(), this->global_blocked_dofs.end()); auto last = std::unique(this->global_blocked_dofs.begin(), this->global_blocked_dofs.end()); this->global_blocked_dofs.resize(last - this->global_blocked_dofs.begin()); auto are_equal = global_blocked_dofs.size() == previous_global_blocked_dofs.size() and std::equal(global_blocked_dofs.begin(), global_blocked_dofs.end(), previous_global_blocked_dofs.begin()); if (not are_equal) { ++this->global_blocked_dofs_release; } } /* -------------------------------------------------------------------------- */ void DOFManager::applyBoundary(const ID & matrix_id) { auto & J = this->getMatrix(matrix_id); if (this->jacobian_release == J.getRelease()) { auto are_equal = this->global_blocked_dofs_release == this->previous_global_blocked_dofs_release; // std::equal(global_blocked_dofs.begin(), global_blocked_dofs.end(), // previous_global_blocked_dofs.begin()); if (not are_equal) { J.applyBoundary(); } previous_global_blocked_dofs.copy(global_blocked_dofs); } else { J.applyBoundary(); } this->jacobian_release = J.getRelease(); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleMatMulVectToGlobalArray(const ID & dof_id, const ID & A_id, const Array & x, SolverVector & array, Real scale_factor) { auto & A = this->getMatrix(A_id); data_cache->zero(); this->assembleToGlobalArray(dof_id, x, *data_cache, 1.); A.matVecMul(*data_cache, array, scale_factor, 1.); } /* -------------------------------------------------------------------------- */ void DOFManager::assembleMatMulVectToResidual(const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor) { assembleMatMulVectToGlobalArray(dof_id, A_id, x, *residual, scale_factor); } } // namespace akantu diff --git a/src/model/common/dof_manager/dof_manager.hh b/src/model/common/dof_manager/dof_manager.hh index 1dfd7a573..08060b2b4 100644 --- a/src/model/common/dof_manager/dof_manager.hh +++ b/src/model/common/dof_manager/dof_manager.hh @@ -1,715 +1,718 @@ /** * @file dof_manager.hh * * @author Nicolas Richart * * @date creation: Tue Aug 18 2015 * @date last modification: Wed Feb 21 2018 * * @brief Class handling the different types of dofs * * * Copyright (©) 2015-2018 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_factory.hh" #include "aka_memory.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #ifndef AKANTU_DOF_MANAGER_HH_ #define AKANTU_DOF_MANAGER_HH_ namespace akantu { class TermsToAssemble; class NonLinearSolver; class TimeStepSolver; class SparseMatrix; class SolverVector; class SolverCallback; } // namespace akantu namespace akantu { class DOFManager : protected Memory, protected MeshEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ protected: struct DOFData; public: DOFManager(const ID & id = "dof_manager", const MemoryID & memory_id = 0); DOFManager(Mesh & mesh, const ID & id = "dof_manager", const MemoryID & memory_id = 0); ~DOFManager() override; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// register an array of degree of freedom virtual void registerDOFs(const ID & dof_id, Array & dofs_array, const DOFSupportType & support_type); /// the dof as an implied type of _dst_nodal and is defined only on a subset /// of nodes virtual void registerDOFs(const ID & dof_id, Array & dofs_array, const ID & support_group); /// register an array of previous values of the degree of freedom virtual void registerDOFsPrevious(const ID & dof_id, Array & dofs_array); /// register an array of increment of degree of freedom virtual void registerDOFsIncrement(const ID & dof_id, Array & dofs_array); /// register an array of derivatives for a particular dof array virtual void registerDOFsDerivative(const ID & dof_id, UInt order, Array & dofs_derivative); /// register array representing the blocked degree of freedoms virtual void registerBlockedDOFs(const ID & dof_id, Array & blocked_dofs); /// Assemble an array to the global residual array virtual void assembleToResidual(const ID & dof_id, Array & array_to_assemble, Real scale_factor = 1.); /// Assemble an array to the global lumped matrix array virtual void assembleToLumpedMatrix(const ID & dof_id, Array & array_to_assemble, const ID & lumped_mtx, Real scale_factor = 1.); /** * Assemble elementary values to a local array of the size nb_nodes * * nb_dof_per_node. The dof number is implicitly considered as * conn(el, n) * nb_nodes_per_element + d. * With 0 < n < nb_nodes_per_element and 0 < d < nb_dof_per_node **/ virtual void assembleElementalArrayLocalArray( const Array & elementary_vect, Array & array_assembeled, ElementType type, GhostType ghost_type, Real scale_factor = 1., const Array & filter_elements = empty_filter); /** * Assemble elementary values to the global residual array. The dof number is * implicitly considered as conn(el, n) * nb_nodes_per_element + d. * With 0 < n < nb_nodes_per_element and 0 < d < nb_dof_per_node **/ virtual void assembleElementalArrayToResidual( const ID & dof_id, const Array & elementary_vect, ElementType type, GhostType ghost_type, Real scale_factor = 1., const Array & filter_elements = empty_filter); /** * Assemble elementary values to a global array corresponding to a lumped * matrix */ virtual void assembleElementalArrayToLumpedMatrix( const ID & dof_id, const Array & elementary_vect, const ID & lumped_mtx, ElementType type, GhostType ghost_type, Real scale_factor = 1., const Array & filter_elements = empty_filter); /** * Assemble elementary values to the global residual array. The dof number is * implicitly considered as conn(el, n) * nb_nodes_per_element + d. With 0 < * n < nb_nodes_per_element and 0 < d < nb_dof_per_node **/ virtual void assembleElementalMatricesToMatrix( const ID & matrix_id, const ID & dof_id, const Array & elementary_mat, ElementType type, GhostType ghost_type = _not_ghost, const MatrixType & elemental_matrix_type = _symmetric, const Array & filter_elements = empty_filter) = 0; /// multiply a vector by a matrix and assemble the result to the residual virtual void assembleMatMulVectToArray(const ID & dof_id, const ID & A_id, const Array & x, Array & array, Real scale_factor = 1) = 0; /// multiply a vector by a lumped matrix and assemble the result to the /// residual virtual void assembleLumpedMatMulVectToResidual(const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor = 1) = 0; /// assemble coupling terms between to dofs virtual void assemblePreassembledMatrix(const ID & dof_id_m, const ID & dof_id_n, const ID & matrix_id, const TermsToAssemble & terms) = 0; /// multiply a vector by a matrix and assemble the result to the residual virtual void assembleMatMulVectToResidual(const ID & dof_id, const ID & A_id, const Array & x, Real scale_factor = 1); /// multiply the dofs by a matrix and assemble the result to the residual virtual void assembleMatMulDOFsToResidual(const ID & A_id, Real scale_factor = 1); /// updates the global blocked_dofs array virtual void updateGlobalBlockedDofs(); /// sets the residual to 0 virtual void zeroResidual(); /// sets the matrix to 0 virtual void zeroMatrix(const ID & mtx); /// sets the lumped matrix to 0 virtual void zeroLumpedMatrix(const ID & mtx); virtual void applyBoundary(const ID & matrix_id = "J"); // virtual void applyBoundaryLumped(const ID & matrix_id = "J"); /// extract a lumped matrix part corresponding to a given dof virtual void getLumpedMatrixPerDOFs(const ID & dof_id, const ID & lumped_mtx, Array & lumped); /// splits the solution storage from a global view to the per dof storages void splitSolutionPerDOFs(); private: /// dispatch the creation of the dof data and register it DOFData & getNewDOFDataInternal(const ID & dof_id); protected: /// common function to help registering dofs the return values are the add new /// numbers of local dofs, pure local dofs, and system size virtual std::tuple registerDOFsInternal(const ID & dof_id, Array & dofs_array); /// minimum functionality to implement per derived version of the DOFManager /// to allow the splitSolutionPerDOFs function to work virtual void getSolutionPerDOFs(const ID & dof_id, Array & solution_array); /// fill a Vector with the equation numbers corresponding to the given /// connectivity static inline void extractElementEquationNumber( const Array & equation_numbers, const Vector & connectivity, UInt nb_degree_of_freedom, Vector & element_equation_number); /// Assemble a array to a global one void assembleMatMulVectToGlobalArray(const ID & dof_id, const ID & A_id, const Array & x, SolverVector & array, Real scale_factor = 1.); /// common function that can be called by derived class with proper matrice /// types template void assemblePreassembledMatrix_(Mat & A, const ID & dof_id_m, const ID & dof_id_n, const TermsToAssemble & terms); template void assembleElementalMatricesToMatrix_( Mat & A, const ID & dof_id, const Array & elementary_mat, ElementType type, GhostType ghost_type, const MatrixType & elemental_matrix_type, const Array & filter_elements); template void assembleMatMulVectToArray_(const ID & dof_id, const ID & A_id, const Array & x, Array & array, Real scale_factor); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// Get the location type of a given dof inline bool isLocalOrMasterDOF(UInt local_dof_num); + /// Is dof periodic (master or slave) + inline bool isPeriodicDOF(UInt local_dof_num); + /// Answer to the question is a dof a slave dof ? inline bool isSlaveDOF(UInt local_dof_num); /// Answer to the question is a dof a slave dof ? inline bool isPureGhostDOF(UInt local_dof_num); /// tells if the dof manager knows about a global dof bool hasGlobalEquationNumber(Int global) const; /// return the local index of the global equation number inline Int globalToLocalEquationNumber(Int global) const; /// converts local equation numbers to global equation numbers; inline Int localToGlobalEquationNumber(Int local) const; /// get the array of dof types (use only if you know what you do...) inline NodeFlag getDOFFlag(Int local_id) const; /// Global number of dofs AKANTU_GET_MACRO(SystemSize, this->system_size, UInt); /// Local number of dofs AKANTU_GET_MACRO(LocalSystemSize, this->local_system_size, UInt); /// Pure local number of dofs AKANTU_GET_MACRO(PureLocalSystemSize, this->pure_local_system_size, UInt); /// Retrieve all the registered DOFs std::vector getDOFIDs() const; /* ------------------------------------------------------------------------ */ /* DOFs and derivatives accessors */ /* ------------------------------------------------------------------------ */ /// Get a reference to the registered dof array for a given id inline Array & getDOFs(const ID & dofs_id); /// Get the support type of a given dof inline DOFSupportType getSupportType(const ID & dofs_id) const; /// are the dofs registered inline bool hasDOFs(const ID & dof_id) const; /// Get a reference to the registered dof derivatives array for a given id inline Array & getDOFsDerivatives(const ID & dofs_id, UInt order); /// Does the dof has derivatives inline bool hasDOFsDerivatives(const ID & dofs_id, UInt order) const; /// Get a reference to the blocked dofs array registered for the given id inline const Array & getBlockedDOFs(const ID & dofs_id) const; /// Does the dof has a blocked array inline bool hasBlockedDOFs(const ID & dofs_id) const; /// Get a reference to the registered dof increment array for a given id inline Array & getDOFsIncrement(const ID & dofs_id); /// Does the dof has a increment array inline bool hasDOFsIncrement(const ID & dofs_id) const; /// Does the dof has a previous array inline Array & getPreviousDOFs(const ID & dofs_id); /// Get a reference to the registered dof array for previous step values a /// given id inline bool hasPreviousDOFs(const ID & dofs_id) const; /// saves the values from dofs to previous dofs virtual void savePreviousDOFs(const ID & dofs_id); /// Get a reference to the solution array registered for the given id inline const Array & getSolution(const ID & dofs_id) const; /// Get a reference to the solution array registered for the given id inline Array & getSolution(const ID & dofs_id); /// Get the blocked dofs array AKANTU_GET_MACRO(GlobalBlockedDOFs, global_blocked_dofs, const Array &); /// Get the blocked dofs array AKANTU_GET_MACRO(PreviousGlobalBlockedDOFs, previous_global_blocked_dofs, const Array &); /* ------------------------------------------------------------------------ */ /* Matrices accessors */ /* ------------------------------------------------------------------------ */ /// Get an instance of a new SparseMatrix virtual SparseMatrix & getNewMatrix(const ID & matrix_id, const MatrixType & matrix_type) = 0; /// Get an instance of a new SparseMatrix as a copy of the SparseMatrix /// matrix_to_copy_id virtual SparseMatrix & getNewMatrix(const ID & matrix_id, const ID & matrix_to_copy_id) = 0; /// Get the equation numbers corresponding to a dof_id. This might be used to /// access the matrix. inline const Array & getLocalEquationsNumbers(const ID & dof_id) const; protected: /// get the array of dof types (use only if you know what you do...) inline const Array & getDOFsAssociatedNodes(const ID & dof_id) const; protected: /* ------------------------------------------------------------------------ */ /// register a matrix SparseMatrix & registerSparseMatrix(const ID & matrix_id, std::unique_ptr & matrix); /// register a lumped matrix (aka a Vector) SolverVector & registerLumpedMatrix(const ID & matrix_id, std::unique_ptr & matrix); /// register a non linear solver instantiated by a derived class NonLinearSolver & registerNonLinearSolver(const ID & non_linear_solver_id, std::unique_ptr & non_linear_solver); /// register a time step solver instantiated by a derived class TimeStepSolver & registerTimeStepSolver(const ID & time_step_solver_id, std::unique_ptr & time_step_solver); template NonLinearSolver & registerNonLinearSolver(DMType & dm, const ID & id, const NonLinearSolverType & type) { ID non_linear_solver_id = this->id + ":nls:" + id; std::unique_ptr nls = std::make_unique( dm, type, non_linear_solver_id, this->memory_id); return this->registerNonLinearSolver(non_linear_solver_id, nls); } template TimeStepSolver & registerTimeStepSolver(DMType & dm, const ID & id, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver, SolverCallback & solver_callback) { ID time_step_solver_id = this->id + ":tss:" + id; std::unique_ptr tss = std::make_unique(dm, type, non_linear_solver, solver_callback, time_step_solver_id, this->memory_id); return this->registerTimeStepSolver(time_step_solver_id, tss); } template SparseMatrix & registerSparseMatrix(DMType & dm, const ID & id, const MatrixType & matrix_type) { ID matrix_id = this->id + ":mtx:" + id; std::unique_ptr sm = std::make_unique(dm, matrix_type, matrix_id); return this->registerSparseMatrix(matrix_id, sm); } template SparseMatrix & registerSparseMatrix(const ID & id, const ID & matrix_to_copy_id) { ID matrix_id = this->id + ":mtx:" + id; auto & sm_to_copy = aka::as_type(this->getMatrix(matrix_to_copy_id)); std::unique_ptr sm = std::make_unique(sm_to_copy, matrix_id); return this->registerSparseMatrix(matrix_id, sm); } template SolverVector & registerLumpedMatrix(DMType & dm, const ID & id) { ID matrix_id = this->id + ":lumped_mtx:" + id; std::unique_ptr sm = std::make_unique(dm, matrix_id); return this->registerLumpedMatrix(matrix_id, sm); } protected: virtual void makeConsistentForPeriodicity(const ID & dof_id, SolverVector & array) = 0; virtual void assembleToGlobalArray(const ID & dof_id, const Array & array_to_assemble, SolverVector & global_array, Real scale_factor) = 0; public: /// extract degrees of freedom (identified by ID) from a global solver array virtual void getArrayPerDOFs(const ID & dof_id, const SolverVector & global, Array & local) = 0; /// Get the reference of an existing matrix SparseMatrix & getMatrix(const ID & matrix_id); /// check if the given matrix exists bool hasMatrix(const ID & matrix_id) const; /// Get an instance of a new lumped matrix virtual SolverVector & getNewLumpedMatrix(const ID & matrix_id) = 0; /// Get the lumped version of a given matrix const SolverVector & getLumpedMatrix(const ID & matrix_id) const; /// Get the lumped version of a given matrix SolverVector & getLumpedMatrix(const ID & matrix_id); /// check if the given matrix exists bool hasLumpedMatrix(const ID & matrix_id) const; /* ------------------------------------------------------------------------ */ /* Non linear system solver */ /* ------------------------------------------------------------------------ */ /// Get instance of a non linear solver virtual NonLinearSolver & getNewNonLinearSolver( const ID & nls_solver_id, const NonLinearSolverType & _non_linear_solver_type) = 0; /// get instance of a non linear solver virtual NonLinearSolver & getNonLinearSolver(const ID & nls_solver_id); /// check if the given solver exists bool hasNonLinearSolver(const ID & solver_id) const; /* ------------------------------------------------------------------------ */ /* Time-Step Solver */ /* ------------------------------------------------------------------------ */ /// Get instance of a time step solver virtual TimeStepSolver & getNewTimeStepSolver(const ID & time_step_solver_id, const TimeStepSolverType & type, NonLinearSolver & non_linear_solver, SolverCallback & solver_callback) = 0; /// get instance of a time step solver virtual TimeStepSolver & getTimeStepSolver(const ID & time_step_solver_id); /// check if the given solver exists bool hasTimeStepSolver(const ID & solver_id) const; /* ------------------------------------------------------------------------ */ const Mesh & getMesh() { if (mesh != nullptr) { return *mesh; } AKANTU_EXCEPTION("No mesh registered in this dof manager"); } /* ------------------------------------------------------------------------ */ AKANTU_GET_MACRO(Communicator, communicator, const auto &); AKANTU_GET_MACRO_NOT_CONST(Communicator, communicator, auto &); /* ------------------------------------------------------------------------ */ AKANTU_GET_MACRO(Solution, *(solution.get()), const auto &); AKANTU_GET_MACRO_NOT_CONST(Solution, *(solution.get()), auto &); AKANTU_GET_MACRO(Residual, *(residual.get()), const auto &); AKANTU_GET_MACRO_NOT_CONST(Residual, *(residual.get()), auto &); /* ------------------------------------------------------------------------ */ /* MeshEventHandler interface */ /* ------------------------------------------------------------------------ */ protected: friend class GlobalDOFInfoDataAccessor; /// helper function for the DOFManager::onNodesAdded method virtual std::pair updateNodalDOFs(const ID & dof_id, const Array & nodes_list); template auto countDOFsForNodes(const DOFData & dof_data, UInt nb_nodes, Func && getNode); void updateDOFsData(DOFData & dof_data, UInt nb_new_local_dofs, UInt nb_new_pure_local, UInt nb_nodes, const std::function & getNode); void updateDOFsData(DOFData & dof_data, UInt nb_new_local_dofs, UInt nb_new_pure_local); auto computeFirstDOFIDs(UInt nb_new_local_dofs, UInt nb_new_pure_local); /// resize all the global information and takes the needed measure like /// cleaning matrices profiles virtual void resizeGlobalArrays(); public: /// function to implement to react on akantu::NewNodesEvent void onNodesAdded(const Array & nodes_list, const NewNodesEvent & event) override; /// function to implement to react on akantu::RemovedNodesEvent void onNodesRemoved(const Array & nodes_list, const Array & new_numbering, const RemovedNodesEvent & event) override; /// function to implement to react on akantu::NewElementsEvent void onElementsAdded(const Array & elements_list, const NewElementsEvent & event) override; /// function to implement to react on akantu::RemovedElementsEvent void onElementsRemoved(const Array & elements_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) override; /// function to implement to react on akantu::ChangedElementsEvent void onElementsChanged(const Array & old_elements_list, const Array & new_elements_list, const ElementTypeMapArray & new_numbering, const ChangedElementsEvent & event) override; protected: inline DOFData & getDOFData(const ID & dof_id); inline const DOFData & getDOFData(const ID & dof_id) const; template inline DOFData_ & getDOFDataTyped(const ID & dof_id); template inline const DOFData_ & getDOFDataTyped(const ID & dof_id) const; virtual std::unique_ptr getNewDOFData(const ID & dof_id) = 0; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// dof representations in the dof manager struct DOFData { DOFData() = delete; explicit DOFData(const ID & dof_id); virtual ~DOFData(); /// DOF support type (nodal, general) this is needed to determine how the /// dof are shared among processors DOFSupportType support_type; ID group_support; /// Degree of freedom array Array * dof{nullptr}; /// Blocked degree of freedoms array Array * blocked_dofs{nullptr}; /// Degree of freedoms increment Array * increment{nullptr}; /// Degree of freedoms at previous step Array * previous{nullptr}; /// Solution associated to the dof Array solution; /* ---------------------------------------------------------------------- */ /* data for dynamic simulations */ /* ---------------------------------------------------------------------- */ /// Degree of freedom derivatives arrays std::vector *> dof_derivatives; /* ---------------------------------------------------------------------- */ /// number of dofs to consider locally for this dof id UInt local_nb_dofs{0}; /// Number of purely local dofs UInt pure_local_nb_dofs{0}; /// number of ghost dofs UInt ghosts_nb_dofs{0}; /// local numbering equation numbers Array local_equation_number; /// associated node for _dst_nodal dofs only Array associated_nodes; virtual Array & getLocalEquationsNumbers() { return local_equation_number; } }; /// type to store dofs information using DOFStorage = std::map>; /// type to store all the matrices using SparseMatricesMap = std::map>; /// type to store all the lumped matrices using LumpedMatricesMap = std::map>; /// type to store all the non linear solver using NonLinearSolversMap = std::map>; /// type to store all the time step solver using TimeStepSolversMap = std::map>; /// store a reference to the dof arrays DOFStorage dofs; /// list of sparse matrices that where created SparseMatricesMap matrices; /// list of lumped matrices LumpedMatricesMap lumped_matrices; /// non linear solvers storage NonLinearSolversMap non_linear_solvers; /// time step solvers storage TimeStepSolversMap time_step_solvers; /// reference to the underlying mesh Mesh * mesh{nullptr}; /// Total number of degrees of freedom (size with the ghosts) UInt local_system_size{0}; /// Number of purely local dofs (size without the ghosts) UInt pure_local_system_size{0}; /// Total number of degrees of freedom UInt system_size{0}; /// rhs to the system of equation corresponding to the residual linked to the /// different dofs std::unique_ptr residual; /// solution of the system of equation corresponding to the different dofs std::unique_ptr solution; /// a vector that helps internally to perform some tasks std::unique_ptr data_cache; /// define the dofs type, local, shared, ghost Array dofs_flag; /// equation number in global numbering Array global_equation_number; using equation_numbers_map = std::unordered_map; /// dual information of global_equation_number equation_numbers_map global_to_local_mapping; /// Communicator used for this manager, should be the same as in the mesh if a /// mesh is registered Communicator & communicator; /// accumulator to know what would be the next global id to use UInt first_global_dof_id{0}; /// Release at last apply boundary on jacobian UInt jacobian_release{0}; /// blocked degree of freedom in the system equation corresponding to the /// different dofs Array global_blocked_dofs; UInt global_blocked_dofs_release{0}; /// blocked degree of freedom in the system equation corresponding to the /// different dofs Array previous_global_blocked_dofs; UInt previous_global_blocked_dofs_release{0}; private: /// This is for unit testing friend class DOFManagerTester; }; using DefaultDOFManagerFactory = Factory; using DOFManagerFactory = Factory; } // namespace akantu #include "dof_manager_inline_impl.hh" #endif /* AKANTU_DOF_MANAGER_HH_ */ diff --git a/src/model/common/dof_manager/dof_manager_inline_impl.hh b/src/model/common/dof_manager/dof_manager_inline_impl.hh index 22713a58d..79b47029c 100644 --- a/src/model/common/dof_manager/dof_manager_inline_impl.hh +++ b/src/model/common/dof_manager/dof_manager_inline_impl.hh @@ -1,335 +1,341 @@ /** * @file dof_manager_inline_impl.hh * * @author Nicolas Richart * * @date creation: Thu Feb 21 2013 * @date last modification: Wed Jan 31 2018 * * @brief inline functions of the dof manager * * * Copyright (©) 2014-2018 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 "dof_manager.hh" #include "element_group.hh" #include "solver_vector.hh" #include "sparse_matrix.hh" #include "terms_to_assemble.hh" /* -------------------------------------------------------------------------- */ #ifndef AKANTU_DOF_MANAGER_INLINE_IMPL_HH_ #define AKANTU_DOF_MANAGER_INLINE_IMPL_HH_ namespace akantu { /* -------------------------------------------------------------------------- */ inline bool DOFManager::hasDOFs(const ID & dof_id) const { auto it = this->dofs.find(dof_id); return it != this->dofs.end(); } /* -------------------------------------------------------------------------- */ inline DOFManager::DOFData & DOFManager::getDOFData(const ID & dof_id) { auto it = this->dofs.find(dof_id); if (it == this->dofs.end()) { AKANTU_EXCEPTION("The dof " << dof_id << " does not exists in " << this->id); } return *it->second; } /* -------------------------------------------------------------------------- */ const DOFManager::DOFData & DOFManager::getDOFData(const ID & dof_id) const { auto it = this->dofs.find(dof_id); if (it == this->dofs.end()) { AKANTU_EXCEPTION("The dof " << dof_id << " does not exists in " << this->id); } return *it->second; } /* -------------------------------------------------------------------------- */ inline void DOFManager::extractElementEquationNumber( const Array & equation_numbers, const Vector & connectivity, UInt nb_degree_of_freedom, Vector & element_equation_number) { for (UInt i = 0, ld = 0; i < connectivity.size(); ++i) { UInt n = connectivity(i); for (UInt d = 0; d < nb_degree_of_freedom; ++d, ++ld) { element_equation_number(ld) = equation_numbers(n * nb_degree_of_freedom + d); } } } /* -------------------------------------------------------------------------- */ template inline DOFData_ & DOFManager::getDOFDataTyped(const ID & dof_id) { return aka::as_type(this->getDOFData(dof_id)); } /* -------------------------------------------------------------------------- */ template inline const DOFData_ & DOFManager::getDOFDataTyped(const ID & dof_id) const { return aka::as_type(this->getDOFData(dof_id)); } /* -------------------------------------------------------------------------- */ inline Array & DOFManager::getDOFs(const ID & dofs_id) { return *(this->getDOFData(dofs_id).dof); } /* -------------------------------------------------------------------------- */ inline DOFSupportType DOFManager::getSupportType(const ID & dofs_id) const { return this->getDOFData(dofs_id).support_type; } /* -------------------------------------------------------------------------- */ inline Array & DOFManager::getPreviousDOFs(const ID & dofs_id) { return *(this->getDOFData(dofs_id).previous); } /* -------------------------------------------------------------------------- */ inline bool DOFManager::hasPreviousDOFs(const ID & dofs_id) const { return (this->getDOFData(dofs_id).previous != nullptr); } /* -------------------------------------------------------------------------- */ inline Array & DOFManager::getDOFsIncrement(const ID & dofs_id) { return *(this->getDOFData(dofs_id).increment); } /* -------------------------------------------------------------------------- */ inline bool DOFManager::hasDOFsIncrement(const ID & dofs_id) const { return (this->getDOFData(dofs_id).increment != nullptr); } /* -------------------------------------------------------------------------- */ inline Array & DOFManager::getDOFsDerivatives(const ID & dofs_id, UInt order) { if (order == 0) { return getDOFs(dofs_id); } std::vector *> & derivatives = this->getDOFData(dofs_id).dof_derivatives; if ((order > derivatives.size()) || (derivatives[order - 1] == nullptr)) { AKANTU_EXCEPTION("No derivatives of order " << order << " present in " << this->id << " for dof " << dofs_id); } return *derivatives[order - 1]; } /* -------------------------------------------------------------------------- */ inline bool DOFManager::hasDOFsDerivatives(const ID & dofs_id, UInt order) const { const std::vector *> & derivatives = this->getDOFData(dofs_id).dof_derivatives; return ((order < derivatives.size()) && (derivatives[order - 1] != nullptr)); } /* -------------------------------------------------------------------------- */ inline const Array & DOFManager::getSolution(const ID & dofs_id) const { return this->getDOFData(dofs_id).solution; } /* -------------------------------------------------------------------------- */ inline Array & DOFManager::getSolution(const ID & dofs_id) { return this->getDOFData(dofs_id).solution; } /* -------------------------------------------------------------------------- */ inline const Array & DOFManager::getBlockedDOFs(const ID & dofs_id) const { return *(this->getDOFData(dofs_id).blocked_dofs); } /* -------------------------------------------------------------------------- */ inline bool DOFManager::hasBlockedDOFs(const ID & dofs_id) const { return (this->getDOFData(dofs_id).blocked_dofs != nullptr); } /* -------------------------------------------------------------------------- */ inline bool DOFManager::isLocalOrMasterDOF(UInt dof_num) { auto dof_flag = this->dofs_flag(dof_num); return (dof_flag & NodeFlag::_local_master_mask) == NodeFlag::_normal; } +/* -------------------------------------------------------------------------- */ +inline bool DOFManager::isPeriodicDOF(UInt dof_num) { + auto dof_flag = this->dofs_flag(dof_num); + return (dof_flag & NodeFlag::_periodic_mask) == NodeFlag::_periodic; +} + /* -------------------------------------------------------------------------- */ inline bool DOFManager::isSlaveDOF(UInt dof_num) { auto dof_flag = this->dofs_flag(dof_num); return (dof_flag & NodeFlag::_shared_mask) == NodeFlag::_slave; } /* -------------------------------------------------------------------------- */ inline bool DOFManager::isPureGhostDOF(UInt dof_num) { auto dof_flag = this->dofs_flag(dof_num); return (dof_flag & NodeFlag::_shared_mask) == NodeFlag::_pure_ghost; } /* -------------------------------------------------------------------------- */ inline Int DOFManager::localToGlobalEquationNumber(Int local) const { return this->global_equation_number(local); } /* -------------------------------------------------------------------------- */ inline bool DOFManager::hasGlobalEquationNumber(Int global) const { auto it = this->global_to_local_mapping.find(global); return (it != this->global_to_local_mapping.end()); } /* -------------------------------------------------------------------------- */ inline Int DOFManager::globalToLocalEquationNumber(Int global) const { auto it = this->global_to_local_mapping.find(global); AKANTU_DEBUG_ASSERT(it != this->global_to_local_mapping.end(), "This global equation number " << global << " does not exists in " << this->id); return it->second; } /* -------------------------------------------------------------------------- */ inline NodeFlag DOFManager::getDOFFlag(Int local_id) const { return this->dofs_flag(local_id); } /* -------------------------------------------------------------------------- */ inline const Array & DOFManager::getDOFsAssociatedNodes(const ID & dof_id) const { const auto & dof_data = this->getDOFData(dof_id); return dof_data.associated_nodes; } /* -------------------------------------------------------------------------- */ const Array & DOFManager::getLocalEquationsNumbers(const ID & dof_id) const { return getDOFData(dof_id).local_equation_number; } /* -------------------------------------------------------------------------- */ template void DOFManager::assembleMatMulVectToArray_(const ID & dof_id, const ID & A_id, const Array & x, Array & array, Real scale_factor) { Vec tmp_array(aka::as_type(*data_cache), this->id + ":tmp_array"); tmp_array.zero(); assembleMatMulVectToGlobalArray(dof_id, A_id, x, tmp_array, scale_factor); getArrayPerDOFs(dof_id, tmp_array, array); } /* -------------------------------------------------------------------------- */ template void DOFManager::assembleElementalMatricesToMatrix_( Mat & A, const ID & dof_id, const Array & elementary_mat, ElementType type, GhostType ghost_type, const MatrixType & elemental_matrix_type, const Array & filter_elements) { AKANTU_DEBUG_IN(); auto & dof_data = this->getDOFData(dof_id); AKANTU_DEBUG_ASSERT(dof_data.support_type == _dst_nodal, "This function applies only on Nodal dofs"); const auto & equation_number = this->getLocalEquationsNumbers(dof_id); UInt nb_element; UInt * filter_it = nullptr; if (filter_elements != empty_filter) { nb_element = filter_elements.size(); filter_it = filter_elements.storage(); } else { if (dof_data.group_support != "__mesh__") { const auto & group_elements = this->mesh->getElementGroup(dof_data.group_support) .getElements(type, ghost_type); nb_element = group_elements.size(); filter_it = group_elements.storage(); } else { nb_element = this->mesh->getNbElement(type, ghost_type); } } AKANTU_DEBUG_ASSERT(elementary_mat.size() == nb_element, "The vector elementary_mat(" << elementary_mat.getID() << ") has not the good size."); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_degree_of_freedom = dof_data.dof->getNbComponent(); const Array & connectivity = this->mesh->getConnectivity(type, ghost_type); auto conn_begin = connectivity.begin(nb_nodes_per_element); auto conn_it = conn_begin; auto size_mat = nb_nodes_per_element * nb_degree_of_freedom; Vector element_eq_nb(nb_degree_of_freedom * nb_nodes_per_element); auto el_mat_it = elementary_mat.begin(size_mat, size_mat); for (UInt e = 0; e < nb_element; ++e, ++el_mat_it) { if (filter_it) { conn_it = conn_begin + *filter_it; } this->extractElementEquationNumber(equation_number, *conn_it, nb_degree_of_freedom, element_eq_nb); std::transform(element_eq_nb.begin(), element_eq_nb.end(), element_eq_nb.begin(), [&](auto && local) { return this->localToGlobalEquationNumber(local); }); if (filter_it) { ++filter_it; } else { ++conn_it; } A.addValues(element_eq_nb, element_eq_nb, *el_mat_it, elemental_matrix_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void DOFManager::assemblePreassembledMatrix_(Mat & A, const ID & dof_id_m, const ID & dof_id_n, const TermsToAssemble & terms) { const auto & equation_number_m = this->getLocalEquationsNumbers(dof_id_m); const auto & equation_number_n = this->getLocalEquationsNumbers(dof_id_n); for (const auto & term : terms) { auto gi = this->localToGlobalEquationNumber(equation_number_m(term.i())); auto gj = this->localToGlobalEquationNumber(equation_number_n(term.j())); A.add(gi, gj, term); } } /* -------------------------------------------------------------------------- */ } // namespace akantu #endif /* AKANTU_DOF_MANAGER_INLINE_IMPL_HH_ */