DomainDecomposition.h

Go to the documentation of this file.

// $Id$
//==============================================================================
//==============================================================================
 
#ifndef _DOMAIN_DECOMPOSITION_H
#define _DOMAIN_DECOMPOSITION_H
 
#include "Interface.h"
#include <map>
#include <set>
#include <string>
#include <vector>
#include <cstddef>
 
class ASMbase;
class ProcessAdm;
class SAMpatch;
class SIMbase;
 
using IntVec = std::vector<int>; 
using IntMat = std::vector<IntVec>; 
 
 
class DomainDecomposition
{
public:
  class OrientIterator {
  public:
    OrientIterator(const ASMbase* pch, int orient, int lIdx, int basis, int dim = 2);
 
    OrientIterator(const ASMbase* pch, int orient, int lIdx);
 
    IntVec::const_iterator begin() { return nodes.begin(); }
    IntVec::const_iterator end() { return nodes.end(); }
    size_t size() const { return nodes.size(); }
 
  protected:
    IntVec nodes;
  };
 
  class SlaveOrder {
    public:
      explicit SlaveOrder(const DomainDecomposition& dd_) : dd(dd_) {}
      SlaveOrder& operator=(const SlaveOrder&) { return *this; }
      bool operator()(const ASM::Interface& A, const ASM::Interface& B) const
      {
        // order by master owner id first
        if (dd.getPatchOwner(A.master) != dd.getPatchOwner(B.master))
          return dd.getPatchOwner(A.master) < dd.getPatchOwner(B.master);
 
        // then by lower slave owner id
        if (dd.getPatchOwner(A.slave) != dd.getPatchOwner(B.slave))
          return dd.getPatchOwner(A.slave) < dd.getPatchOwner(B.slave);
 
        // then by slave id
        if (A.slave != B.slave)
          return A.slave < B.slave;
 
        // then by dim
        if (A.dim != B.dim)
          return A.dim < B.dim;
 
        // then by basis
        if (A.basis != B.basis)
          return A.basis < B.basis;
 
        // finally by index on master
        return A.midx < B.midx;
      }
    protected:
      const DomainDecomposition& dd; 
  };
 
  std::set<ASM::Interface, SlaveOrder> ghostConnections; 
 
  DomainDecomposition() : ghostConnections(SlaveOrder(*this)), blocks(1) {}
 
  bool setup(const ProcessAdm& adm, const SIMbase& sim);
 
  bool setup(const ProcessAdm& adm, const IntMat& neighs, const IntMat& meqn);
 
  void setup(const int neq);
 
  std::vector<std::set<int>>
  getSubdomains(int nx, int ny, int nz, int overlap, size_t block) const;
 
  static IntMat calcSubdomains(size_t nel1, size_t nel2, size_t nel3,
                               size_t g1, size_t g2, size_t g3, size_t overlap);
 
  int getMinEq(size_t idx = 0) const { return blocks[idx].minEq; }
  int getMaxEq(size_t idx = 0) const { return blocks[idx].maxEq; }
  int getNoGlbEqs(size_t idx = 0) const { return blocks[idx].nGlbEqs; }
  int getMinNode() const { return minNode; }
  int getMaxNode() const { return maxNode; }
  int getMinDOF() const { return minDof; }
  int getMaxDOF() const { return maxDof; }
 
  void setPatchOwner(size_t p, size_t owner) { patchOwner[p] = owner; }
 
  int getPatchOwner(size_t p) const;
 
  int getGlobalEq(int lEq, size_t idx=0) const;
 
  const IntVec& getMLGEQ(size_t idx = 0) const
  { return blocks[idx].MLGEQ; }
 
  const std::set<int>& getBlockEqs(size_t idx) const
  { return blocks[idx+1].localEqs; }
 
  const std::map<int,int>& getG2LEQ(size_t idx) const
  { return blocks[idx].G2LEQ; }
 
  const IntVec& getMLGN() const { return MLGN; }
 
  const SAMpatch* getSAM() const { return sam; }
 
  size_t getNoBlocks() const { return blocks.size()-1; }
 
  bool isPartitioned() const { return !myElms.empty(); }
 
  const IntVec& getElms() const { return myElms; }
 
  void setElms(const IntVec& elms, const std::string& save)
  { myElms = elms; savePart = save; }
 
  void clear();
 
private:
 static IntMat calcSubdomains1D(size_t nel1, size_t g1, size_t overlap);
 
  static IntMat calcSubdomains2D(size_t nel1, size_t nel2,
                                 size_t g1, size_t g2, size_t overlap);
 
  static IntMat calcSubdomains3D(size_t nel1, size_t nel2, size_t nel3,
                                 size_t g1, size_t g2, size_t g3,
                                 size_t overlap);
 
#ifdef HAVE_MPI
  IntVec setupEquationNumbers(const SIMbase& sim,
                              int pidx, int lidx,
                              const std::set<int>& cbasis,
                              int dim, int thick, int orient = 0);
 
  void setupNodeNumbers(int basis, IntVec& lNodes,
                        std::set<int>& cbasis,
                        const ASMbase* pch,
                        int dim, int lidx, int thick, int orient = 0);
#endif
 
  bool calcGlobalNodeNumbers(const ProcessAdm& adm, const SIMbase& sim);
 
  bool calcGlobalEqNumbers(const ProcessAdm& adm, const SIMbase& sim);
 
  template<class Callback>
  bool calcGlobalEqNumbersPart(const ProcessAdm& adm, const Callback& cb);
 
  bool sanityCheckCorners(const SIMbase& sim);
 
  bool graphPartition(const ProcessAdm& adm, const IntMat& neigh);
 
  std::map<int,int> patchOwner; 
 
  struct BlockInfo {
    int basis;      
    int components; 
    IntVec MLGEQ; 
    int minEq; 
    int maxEq; 
    int nGlbEqs; 
    std::set<int> localEqs; 
    std::map<int,int> G2LEQ; 
  };
 
  IntVec MLGN; 
  std::vector<BlockInfo> blocks; 
  IntVec myElms; 
  int minDof = 0; 
  int maxDof = 0; 
  int minNode = 0; 
  int maxNode = 0; 
 
  const SAMpatch* sam = nullptr; 
 
  std::string savePart; 
};
 
#endif