SF_partitioning.h

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// openCARP is an open cardiac electrophysiology simulator.
//
// Copyright (C) 2020 openCARP project
//
// This program is licensed under the openCARP Academic Public License (APL)
// v1.0: You can use and redistribute it and/or modify it in non-commercial
// academic environments under the terms of APL as published by the openCARP
// project v1.0, or (at your option) any later version. Commercial use requires
// a commercial license (info@opencarp.org).
//
// This program is distributed without any warranty; see the openCARP APL for
// more details.
//
// You should have received a copy of the openCARP APL along with this program
// and can find it online: http://www.opencarp.org/license
// ----------------------------------------------------------------------------

#ifndef _SF_PARTITIONING_H
#define _SF_PARTITIONING_H

#define KDPART_MPI
#include "kdpart.hpp"
#include "SF_container.h"
#include "SF_vector.h"

namespace SF {

template<class T, class S>
class abstract_partitioner
{
public:
  virtual void operator()(const meshdata<T, S> & mesh, vector<T> & part) = 0;
};

#ifdef WITH_PARMETIS

#include "parmetis.h"

template<class T, class S>
class parmetis_partitioner : public abstract_partitioner<T, S>
{
private:
  // this are the options configurable by the user
  float _unbalance; 
  short _ncommon;   

public:
  parmetis_partitioner() : _unbalance(1.01f), _ncommon(2)
  {}
  parmetis_partitioner(float ub, short nc) : _unbalance(ub), _ncommon(nc)
  {}

  inline void operator()(const meshdata<T, S> & mesh, vector<T> & part)
  {
    MPI_Comm comm = mesh.comm;

    int rank, size;
    MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);

    part.resize(mesh.l_numelem);
    // special treatment of sequential case since parmetis is rather slow
    if(size == 1) {
      part.zero();
      return;
    }

    // communicate the distribution of the elements
    vector<idx_t> elemdist_cnt(size), elemdist_dsp(size+1);
    idx_t numelem = mesh.l_numelem; // we need a type conversion to int
    MPI_Allgather(&numelem, sizeof(idx_t), MPI_BYTE, elemdist_cnt.data(), sizeof(idx_t), MPI_BYTE, comm);
    dsp_from_cnt(elemdist_cnt, elemdist_dsp);

    int    nparts     = size;     // number of partitions
    int    ncommonnod = _ncommon; // number of common nodes defining the neighbourhood between 2 elems
    int    ncon       = 1; // number of constraints to satisfy (we always use 1)
    int    wgtflag    = 0;        // number of weights for graph edges (unused, thus 0)
    idx_t* elemwgt    = NULL;
    int    numflag    = 0;        // indexing is 0 based

    vector<real_t> tpwgts(nparts*ncon); // weights for nodes
    for(size_t i=0; i<tpwgts.size(); i++) tpwgts[i] = (float)1. / (float)(nparts*ncon);

    vector<real_t> ubvec(ncon);
    for (int i = 0; i < ncon; i++ ) ubvec[i] = _unbalance;

    vector<int> opt(4, 0); // options are unused

    // we need to convert the mesh to idx_t
    vector<idx_t> eptr(mesh.dsp.size()), eind(mesh.con.size());
    vec_assign(eptr.data(), mesh.dsp.data(), eptr.size());

    const vector<T> & rnod = mesh.get_numbering(NBR_REF);
    for(size_t i=0; i<mesh.con.size(); i++) eind[i] = rnod[mesh.con[i]];

    vector<idx_t> pm_part(mesh.l_numelem);

    int n_edgecut;
    ParMETIS_V3_PartMeshKway(elemdist_dsp.data(),
                             eptr.data(),
                             eind.data(),
                             elemwgt,
                             &wgtflag,
                             &numflag,
                             &ncon,
                             &ncommonnod,
                             &nparts,
                             tpwgts.data(),
                             ubvec.data(),
                             opt.data(),
                             &n_edgecut,
                             pm_part.data(),
                             &comm);

    vec_assign(part.data(), pm_part.data(), part.size());

  }
};

#endif

template<class T, class S>
class kdtree_partitioner : public abstract_partitioner<T,S>
{
  public:
  void operator() (const meshdata<T,S> & mesh, vector<T> & part_vec)
  {
    int size; MPI_Comm_size(mesh.comm, &size);

    part_vec.resize(mesh.l_numelem);
    // special treatment of sequential case
    if(size == 1) {
      part_vec.zero();
      return;
    }

    std::vector<S> ctr(mesh.l_numelem*3);
    for(size_t i=0; i<mesh.l_numelem; i++) {
      kdpart::vec3<S> avrg;
      T dsp = mesh.dsp[i];
      T elemsize = mesh.dsp[i+1] - dsp;

      for(T j=0; j<elemsize; j++) {
        T v = mesh.con[dsp+j];
        avrg.x += mesh.xyz[v*3+0];
        avrg.y += mesh.xyz[v*3+1];
        avrg.z += mesh.xyz[v*3+2];
      }
      avrg.x /= S(elemsize);
      avrg.y /= S(elemsize);
      avrg.z /= S(elemsize);

      ctr[i*3+0] = avrg.x;
      ctr[i*3+1] = avrg.y;
      ctr[i*3+2] = avrg.z;
    }

    std::vector<T> part;
    kdpart::parallel_partitioner<T,S> partitioner;
    partitioner(mesh.comm, ctr, size, part);

    part_vec.assign(part.begin(), part.end());
  }
};

}

#endif