pncmesh.hpp

Go to the documentation of this file.

// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.org.
//
// MFEM 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) version 2.1 dated February 1999.

#ifndef MFEM_PNCMESH
#define MFEM_PNCMESH

#include "../config/config.hpp"

#ifdef MFEM_USE_MPI

#include <map>
#include <set>

#include "ncmesh.hpp"
#include "../general/communication.hpp"

namespace mfem
{

class ParMesh;
class FiniteElementCollection; // for edge orientation handling
class FiniteElementSpace; // for Dof -> VDof conversion


class ParNCMesh : public NCMesh
{
public:
   ParNCMesh(MPI_Comm comm, const NCMesh& ncmesh);

   virtual ~ParNCMesh();

   virtual void Refine(const Array<Refinement> &refinements);

   virtual void LimitNCLevel(int max_nc_level);

   virtual void CheckDerefinementNCLevel(const Table &deref_table,
                                         Array<int> &level_ok, int max_nc_level);

   virtual void Derefine(const Array<int> &derefs);

   void Rebalance();


   // interface for ParFiniteElementSpace

   const NCList& GetSharedVertices()
   {
      if (shared_vertices.Empty()) { BuildSharedVertices(); }
      return shared_vertices;
   }

   const NCList& GetSharedEdges()
   {
      if (edge_list.Empty()) { BuildEdgeList(); }
      return shared_edges;
   }

   const NCList& GetSharedFaces()
   {
      if (face_list.Empty()) { BuildFaceList(); }
      return shared_faces;
   }

   const NCList& GetSharedList(int type)
   {
      switch (type)
      {
         case 0: return GetSharedVertices();
         case 1: return GetSharedEdges();
         default: return GetSharedFaces();
      }
   }

   int GetFaceOrientation(int index) const
   {
      return face_orient[index];
   }

   int GetOwner(int type, int index) const
   {
      switch (type)
      {
         case 0: return vertex_owner[index];
         case 1: return edge_owner[index];
         default: return face_owner[index];
      }
   }

   const int* GetGroup(int type, int index, int &size) const
   {
      const Table* table;
      switch (type)
      {
         case 0: table = &vertex_group; break;
         case 1: table = &edge_group; break;
         default: table = &face_group;
      }
      size = table->RowSize(index);
      return table->GetRow(index);
   }

   bool RankInGroup(int type, int index, int rank) const
   {
      int size;
      const int* group = GetGroup(type, index, size);
      for (int i = 0; i < size; i++)
      {
         if (group[i] == rank) { return true; }
      }
      return false;
   }

   bool IsGhost(int type, int index) const
   {
      switch (type)
      {
         case 0: return index >= NVertices;
         case 1: return index >= NEdges;
         default: return index >= NFaces;
      }
   }

   int ElementRank(int index) const
   { return leaf_elements[index]->rank; }


   // interface for ParMesh

   void GetFaceNeighbors(ParMesh &pmesh);


   // utility

   void SendRebalanceDofs(int old_ndofs, const Table &old_element_dofs,
                          long old_global_offset, FiniteElementSpace* space);

   void RecvRebalanceDofs(Array<int> &elements, Array<long> &dofs);

   const Array<int>& GetRebalanceOldIndex() const { return old_index_or_rank; }

   const Array<int>& GetDerefineOldRanks() const { return old_index_or_rank; }

   template<typename Type>
   void SynchronizeDerefinementData(Array<Type> &elem_data,
                                    const Table &deref_table);

   virtual void GetBoundaryClosure(const Array<int> &bdr_attr_is_ess,
                                   Array<int> &bdr_vertices,
                                   Array<int> &bdr_edges);

   void GetDebugMesh(Mesh &debug_mesh) const;


protected:
   MPI_Comm MyComm;
   int NRanks, MyRank;

   int NVertices, NGhostVertices;
   int NEdges, NGhostEdges;
   int NFaces, NGhostFaces;
   int NElements, NGhostElements;

   // lists of vertices/edges/faces shared by us and at least one more processor
   NCList shared_vertices;
   NCList shared_edges;
   NCList shared_faces;

   // owner processor for each vertex/edge/face
   Array<int> vertex_owner;
   Array<int> edge_owner;
   Array<int> face_owner;

   // list of processors sharing each vertex/edge/face
   Table vertex_group;
   Table edge_group;
   Table face_group;

   Array<char> face_orient; // see CalcFaceOrientations

   Array<char> element_type;

   Array<Element*> ghost_layer; 
   Array<Element*> boundary_layer; 

   virtual void Update();

   int Partition(long index, long total_elements) const
   { return index * NRanks / total_elements; }

   int InitialPartition(int index) const
   { return Partition(index, leaf_elements.Size()); }

   long PartitionFirstIndex(int rank, long total_elements) const
   { return (rank * total_elements + NRanks-1) / NRanks; }

   virtual void UpdateVertices();
   virtual void AssignLeafIndices();

   virtual bool IsGhost(const Element* elem) const
   { return elem->rank != MyRank; }

   virtual int GetNumGhosts() const { return NGhostElements; }

   virtual void OnMeshUpdated(Mesh *mesh);

   virtual void BuildEdgeList();
   virtual void BuildFaceList();

   virtual void ElementSharesEdge(Element* elem, Edge* edge);
   virtual void ElementSharesFace(Element* elem, Face* face);

   void BuildSharedVertices();

   static int get_face_orientation(Face *face, Element* e1, Element* e2,
                                   int local[2] = NULL);
   void CalcFaceOrientations();

   void UpdateLayers();

   Array<Connection> index_rank; // temporary

   void AddMasterSlaveRanks(int nitems, const NCList& list);
   void MakeShared(const Table &groups, const NCList &list, NCList &shared);

   class ElementSet
   {
   public:
      ElementSet(NCMesh *ncmesh = NULL, bool include_ref_types = false)
         : ncmesh(ncmesh), include_ref_types(include_ref_types) {}
      ElementSet(const ElementSet &other);

      void Encode(const Array<Element*> &elements);
      void Dump(std::ostream &os) const;

      void Load(std::istream &is);
      void Decode(Array<Element*> &elements) const;

      void SetNCMesh(NCMesh *ncmesh) { this->ncmesh = ncmesh; }

   protected:
      Array<unsigned char> data; 
      NCMesh* ncmesh;
      bool include_ref_types;

      void EncodeTree(Element* elem);
      void DecodeTree(Element* elem, int &pos, Array<Element*> &elements) const;

      void WriteInt(int value);
      int  GetInt(int pos) const;
      void FlagElements(const Array<Element*> &elements, char flag);
   };

   void EncodeMeshIds(std::ostream &os, Array<MeshId> ids[]);

   void DecodeMeshIds(std::istream &is, Array<MeshId> ids[]);

   bool CheckElementType(Element* elem, int type);

   Array<Element*> tmp_neighbors; // temporary used by ElementNeighborProcessors

   void ElementNeighborProcessors(Element* elem, Array<int> &ranks);

   void NeighborProcessors(Array<int> &neighbors);

   void Prune();

   bool PruneTree(Element* elem);


   template<class ValueType, bool RefTypes, int Tag>
   class ElementValueMessage : public VarMessage<Tag>
   {
   public:
      using VarMessage<Tag>::data;
      std::vector<Element*> elements;
      std::vector<ValueType> values;

      int Size() const { return elements.size(); }
      void Reserve(int size) { elements.reserve(size); values.reserve(size); }

      void Add(Element* elem, ValueType val)
      { elements.push_back(elem); values.push_back(val); }

      void SetNCMesh(ParNCMesh* pncmesh) { this->pncmesh = pncmesh; }

      ElementValueMessage() : pncmesh(NULL) {}

   protected:
      ParNCMesh* pncmesh;

      virtual void Encode();
      virtual void Decode();
   };

   class NeighborRefinementMessage : public ElementValueMessage<char, false, 289>
   {
   public:
      void AddRefinement(Element* elem, char ref_type) { Add(elem, ref_type); }
      typedef std::map<int, NeighborRefinementMessage> Map;
   };

   class NeighborDerefinementMessage : public ElementValueMessage<int, false, 290>
   {
   public:
      void AddDerefinement(Element* elem, int rank) { Add(elem, rank); }
      typedef std::map<int, NeighborDerefinementMessage> Map;
   };

   class NeighborElementRankMessage : public ElementValueMessage<int, false, 156>
   {
   public:
      void AddElementRank(Element* elem, int rank) { Add(elem, rank); }
      typedef std::map<int, NeighborElementRankMessage> Map;
   };

   class RebalanceMessage : public ElementValueMessage<int, true, 157>
   {
   public:
      void AddElementRank(Element* elem, int rank) { Add(elem, rank); }
      typedef std::map<int, RebalanceMessage> Map;
   };

   class RebalanceDofMessage : public VarMessage<158>
   {
   public:
      std::vector<int> elem_ids, dofs;
      long dof_offset;

      void SetElements(const Array<Element*> &elems, NCMesh *ncmesh);
      void SetNCMesh(NCMesh* ncmesh) { eset.SetNCMesh(ncmesh); }

      typedef std::map<int, RebalanceDofMessage> Map;

   protected:
      ElementSet eset;

      virtual void Encode();
      virtual void Decode();
   };

   void RedistributeElements(Array<int> &new_ranks, int target_elements,
                             bool record_comm);

   RebalanceDofMessage::Map send_rebalance_dofs;
   RebalanceDofMessage::Map recv_rebalance_dofs;

   Array<int> old_index_or_rank;

   Array<DenseMatrix*> aux_pm_store;
   void ClearAuxPM();

   static bool compare_ranks(const Element* a, const Element* b);
   static bool compare_ranks_indices(const Element* a, const Element* b);

   friend class ParMesh;
   friend class NeighborDofMessage;
};


class NeighborDofMessage : public VarMessage<135>
{
public:
   void AddDofs(int type, const NCMesh::MeshId &id, const Array<int> &dofs);

   void Init(ParNCMesh* pncmesh, const FiniteElementCollection* fec, int ndofs)
   { this->pncmesh = pncmesh; this->fec = fec; this->ndofs = ndofs; }

   void GetDofs(int type, const NCMesh::MeshId& id,
                Array<int>& dofs, int &ndofs);

   typedef std::map<int, NeighborDofMessage> Map;

protected:
   typedef std::map<NCMesh::MeshId, std::vector<int> > IdToDofs;
   IdToDofs id_dofs[3];

   ParNCMesh* pncmesh;
   const FiniteElementCollection* fec;
   int ndofs;

   virtual void Encode();
   virtual void Decode();

   void ReorderEdgeDofs(const NCMesh::MeshId &id, std::vector<int> &dofs);
};

class NeighborRowRequest: public VarMessage<312>
{
public:
   std::set<int> rows;

   void RequestRow(int row) { rows.insert(row); }
   void RemoveRequest(int row) { rows.erase(row); }

   typedef std::map<int, NeighborRowRequest> Map;

protected:
   virtual void Encode();
   virtual void Decode();
};

class NeighborRowReply: public VarMessage<313>
{
public:
   void AddRow(int row, const Array<int> &cols, const Vector &srow);

   bool HaveRow(int row) const { return rows.find(row) != rows.end(); }
   void GetRow(int row, Array<int> &cols, Vector &srow);

   typedef std::map<int, NeighborRowReply> Map;

protected:
   struct Row { std::vector<int> cols; Vector srow; };
   std::map<int, Row> rows;

   virtual void Encode();
   virtual void Decode();
};


// comparison operator so that MeshId can be used as key in std::map
inline bool operator< (const NCMesh::MeshId &a, const NCMesh::MeshId &b)
{
   return a.index < b.index;
}

} // namespace mfem

#endif // MFEM_USE_MPI

#endif // MFEM_PNCMESH