50#ifndef _ZOLTAN2_APFMESHADAPTER_HPP_
51#define _ZOLTAN2_APFMESHADAPTER_HPP_
56#include <unordered_map>
61#ifndef HAVE_ZOLTAN2_PARMA
67template <
typename User>
73 APFMeshAdapter(
const Comm<int> &comm, apf::Mesh* m,std::string primary,std::string adjacency,
bool needSecondAdj=
false)
75 throw std::runtime_error(
76 "BUILD ERROR: ParMA requested but not compiled into Zoltan2.\n"
77 "Please set CMake flag Trilinos_ENABLE_SCOREC:BOOL=ON.");
83#ifdef HAVE_ZOLTAN2_PARMA
86#include <apfDynamicArray.h>
87#include <apfNumbering.h>
118template <
typename User>
119 class APFMeshAdapter:
public MeshAdapter<User> {
143 APFMeshAdapter(
const Comm<int> &comm, apf::Mesh* m,std::string primary,
144 std::string adjacency,
bool needSecondAdj=
false,
int needs=0);
147 void print(
int me,
int verbosity=0);
148 template <
typename Adapter>
150 const PartitioningSolution<Adapter> &solution)
const{
152 apf::Migration* plan =
new apf::Migration(*out);
153 const part_t* new_part_ids = solution.getPartListView();
158 apf::MeshIterator* itr = (*out)->begin(m_dimension);
159 apf::MeshEntity* ent;
161 while ((ent=(*out)->iterate(itr))) {
162 assert(new_part_ids[i]<PCU_Comm_Peers());
163 plan->send(ent,new_part_ids[i]);
172 apf::MeshIterator* itr = (*out)->begin(m_dimension);
173 apf::MeshEntity* ent;
175 while ((ent=(*out)->iterate(itr))) {
176 std::unordered_map<unsigned int,unsigned int> newOwners;
178 unsigned int max_num = 0;
179 int new_part=PCU_Comm_Self();
180 unsigned int num = in->getDownward(ent,dim,adj);
181 for (
unsigned int j=0;j<num;j++) {
182 gno_t gid = apf::getNumber(gids[dim],apf::Node(adj[j],0));
183 lno_t lid = apf::getNumber(lids[dim],adj[j],0,0);
184 newOwners[new_part_ids[lid]]++;
185 if (newOwners[new_part_ids[lid]]>max_num) {
186 max_num=newOwners[new_part_ids[lid]];
187 new_part = new_part_ids[lid];
191 if (new_part<0||new_part>=PCU_Comm_Peers()) {
192 std::cout<<new_part<<std::endl;
193 throw std::runtime_error(
"Target part is out of bounds\n");
195 plan->send(ent,new_part);
200 (*out)->migrate(plan);
214 int dim = entityZ2toAPF(etype);
215 return dim==m_dimension;
219 int dim = entityZ2toAPF(etype);
220 if (dim<=m_dimension&&dim>=0)
221 return num_local[dim];
227 int dim = entityZ2toAPF(etype);
228 if (dim<=m_dimension&&dim>=0)
229 Ids = gid_mapping[dim];
236 int dim = entityZ2toAPF(etype);
237 if (dim<=m_dimension&&dim>=0)
238 Types = topologies[dim];
244 int dim = entityZ2toAPF(etype);
245 return static_cast<int>(
weights[dim].size());
249 int &stride,
int idx = 0)
const
251 int dim = entityZ2toAPF(etype);
252 typename map_array_t::iterator itr =
weights[dim].find(idx);
254 ws = &(*(itr->second.first));
255 stride = itr->second.second;
266 int &stride,
int coordDim)
const {
267 if (coordDim>=0 && coordDim<3) {
268 int dim = entityZ2toAPF(etype);
269 if (dim<=m_dimension&&dim>=0) {
270 coords = ent_coords[dim]+coordDim;
285 int dim_source = entityZ2toAPF(source);
286 int dim_target = entityZ2toAPF(target);
287 return dim_source<=m_dimension && dim_source>=0 &&
288 dim_target<=m_dimension && dim_target>=0 &&
289 dim_target!=dim_source&&
290 has(dim_source) && has(dim_target);
295 int dim_source = entityZ2toAPF(source);
296 int dim_target = entityZ2toAPF(target);
298 return adj_gids[dim_source][dim_target].size();
305 int dim_source = entityZ2toAPF(source);
306 int dim_target = entityZ2toAPF(target);
308 offsets = adj_offsets[dim_source][dim_target];
309 adjacencyIds = &(adj_gids[dim_source][dim_target][0]);
319#ifndef USE_MESH_ADAPTER
324 int dim_source = entityZ2toAPF(sourcetarget);
325 int dim_target = entityZ2toAPF(through);
326 if (dim_source==1&&dim_target==0)
328 return dim_source<=m_dimension && dim_source>=0 &&
329 dim_target<=m_dimension && dim_target>=0 &&
330 dim_target!=dim_source &&
331 has(dim_source)&&has(dim_target);
337 int dim_source = entityZ2toAPF(sourcetarget);
338 int dim_target = entityZ2toAPF(through);
340 return adj2_gids[dim_source][dim_target].size();
348 int dim_source = entityZ2toAPF(sourcetarget);
349 int dim_target = entityZ2toAPF(through);
351 offsets=adj2_offsets[dim_source][dim_target];
352 adjacencyIds=&(adj2_gids[dim_source][dim_target][0]);
391 bool has(
int dim)
const {
return (entity_needs>>dim)%2;}
394 int entityZ2toAPF(
enum MeshEntityType etype)
const {
return static_cast<int>(etype);}
398 if (ttype==apf::Mesh::VERTEX)
400 else if (ttype==apf::Mesh::EDGE)
402 else if (ttype==apf::Mesh::TRIANGLE)
404 else if (ttype==apf::Mesh::QUAD)
406 else if (ttype==apf::Mesh::TET)
408 else if (ttype==apf::Mesh::HEX)
410 else if (ttype==apf::Mesh::PRISM)
412 else if (ttype==apf::Mesh::PYRAMID)
415 throw "No such APF topology type";
420 void getTagWeight(apf::Mesh* m, apf::MeshTag* tag,apf::MeshEntity* ent,
scalar_t* ws);
429 apf::Numbering** lids;
430 apf::GlobalNumbering** gids;
435 std::vector<gno_t>** adj_gids;
437 std::vector<gno_t>** adj2_gids;
442 typedef std::unordered_map<int, std::pair<ArrayRCP<const scalar_t>,
int> > map_array_t;
451template <
typename User>
455 std::string adjacency,
460 m_dimension = m->getDimension();
463 entity_needs = needs;
469 if (primary==
"element") {
475 if (adjacency==
"element") {
481 if (primary==
"region"&&m_dimension<3)
482 throw std::runtime_error(
"primary type and mesh dimension mismatch");
483 if (adjacency==
"region"&&m_dimension<3)
484 throw std::runtime_error(
"adjacency type and mesh dimension mismatch");
485 this->setEntityTypes(primary,adjacency,adjacency);
488 int dim1 = entityZ2toAPF(this->getPrimaryEntityType());
489 int dim2 = entityZ2toAPF(this->getAdjacencyEntityType());
493 entity_needs|=new_needs;
496 lids =
new apf::Numbering*[m_dimension+1];
497 gids =
new apf::GlobalNumbering*[m_dimension+1];
498 gid_mapping =
new gno_t*[m_dimension+1];
499 std::unordered_map<gno_t,lno_t>* lid_mapping =
new std::unordered_map<gno_t,lno_t>[m_dimension+1];
500 num_local =
new size_t[m_dimension+1];
503 for (
int i=0;i<=m_dimension;i++) {
506 topologies[i] = NULL;
507 gid_mapping[i] = NULL;
511 num_local[i] = m->count(i);
512 long global_count = countOwned(m,i);
513 PCU_Add_Longs(&global_count,1);
518 sprintf(lids_name,
"lids%d",i);
520 sprintf(gids_name,
"ids%d",i);
521 apf::FieldShape*
shape = apf::getConstant(i);
522 lids[i] = apf::createNumbering(m,lids_name,shape,1);
523 apf::Numbering* tmp = apf::numberOwnedDimension(m,gids_name,i);
524 gids[i] = apf::makeGlobal(tmp);
525 apf::synchronize(gids[i]);
526 apf::MeshIterator* itr = m->begin(i);
527 apf::MeshEntity* ent;
529 while ((ent=m->iterate(itr))) {
530 apf::number(lids[i],ent,0,0,num);
531 lid_mapping[i][apf::getNumber(gids[i],apf::Node(ent,0))]=num;
535 assert(num==num_local[i]);
539 gid_mapping[i] =
new gno_t[num_local[i]];
541 apf::DynamicArray<apf::Node> nodes;
544 while((ent=m->iterate(itr))) {
545 gno_t gid = apf::getNumber(gids[i],apf::Node(ent,0));
546 gid_mapping[i][ apf::getNumber(lids[i],ent,0,0)] = gid;
555 adj_gids =
new std::vector<gno_t>*[m_dimension+1];
556 adj_offsets =
new offset_t**[m_dimension+1];
558 adj2_gids =
new std::vector<gno_t>*[m_dimension+1];
559 adj2_offsets =
new offset_t**[m_dimension+1];
565 for (
int i=0;i<=m_dimension;i++) {
570 adj2_offsets[i]=NULL;
574 adj_gids[i] =
new std::vector<gno_t>[m_dimension+1];
575 adj_offsets[i] =
new offset_t*[m_dimension+1];
577 adj2_gids[i] =
new std::vector<gno_t>[m_dimension+1];
578 adj2_offsets[i] =
new offset_t*[m_dimension+1];
580 for (
int j=0;j<=m_dimension;j++) {
583 adj_offsets[i][j]=NULL;
585 adj2_offsets[i][j]=NULL;
590 apf::MeshIterator* itr = m->begin(i);
591 apf::MeshEntity* ent;
593 adj_offsets[i][j] =
new offset_t[num_local[i]+1];
594 adj_offsets[i][j][0] =0;
596 adj2_offsets[i][j] =
new offset_t[num_local[i]+1];
597 adj2_offsets[i][j][0] =0;
604 std::unordered_map<gno_t,apf::MeshEntity*> part_boundary_mapping;
606 while ((ent=m->iterate(itr))) {
607 std::set<gno_t> temp_adjs;
610 m->getAdjacent(ent,j,adj);
611 for (
unsigned int l=0;l<adj.getSize();l++) {
612 adj_gids[i][j].push_back(apf::getNumber(gids[j],apf::Node(adj[l],0)));
616 m->getAdjacent(adj[l],i,adj2);
617 for (
unsigned int o=0;o<adj2.getSize();o++)
618 temp_adjs.insert(apf::getNumber(gids[i],apf::Node(adj2[o],0)));
619 if (i==m_dimension) {
621 m->getResidence(adj[l],res);
623 part_boundary_mapping[apf::getNumber(gids[j],apf::Node(adj[l],0))] = adj[l];
624 for (apf::Parts::iterator it=res.begin();it!=res.end();it++) {
626 send_vals[1]=apf::getNumber(gids[i],apf::Node(ent,0));
627 send_vals[0]=apf::getNumber(gids[j],apf::Node(adj[l],0));
629 PCU_Comm_Pack(*it,send_vals,2*
sizeof(
gno_t));
634 adj_offsets[i][j][k] = adj_gids[i][j].size();
637 if (needSecondAdj && i!=m_dimension) {
639 m->getResidence(ent,res);
640 typename std::set<gno_t>::iterator adj_itr;
641 for (adj_itr=temp_adjs.begin();adj_itr!=temp_adjs.end();adj_itr++) {
642 for (apf::Parts::iterator it=res.begin();it!=res.end();it++) {
644 send_vals[0]=apf::getNumber(gids[i],apf::Node(ent,0));
645 send_vals[1] = *adj_itr;
646 if (send_vals[0]!=send_vals[1])
647 PCU_Comm_Pack(*it,send_vals,2*
sizeof(
gno_t));
656 std::set<gno_t>* adjs2 =
new std::set<gno_t>[num_local[i]];
657 while (PCU_Comm_Receive()) {
659 PCU_Comm_Unpack(adj2,2*
sizeof(
gno_t));
661 if (i==m_dimension) {
662 apf::MeshEntity* through = part_boundary_mapping[adj2[0]];
664 m->getAdjacent(through,i,adj);
665 for (
unsigned int l=0;l<adj.getSize();l++) {
666 if (apf::getNumber(gids[i],apf::Node(adj[l],0))!=adj2[1])
667 adjs2[apf::getNumber(lids[i],adj[l],0,0)].insert(adj2[1]);
671 lno_t index = lid_mapping[i][adj2[0]];
672 adjs2[index].insert(adj2[1]);
676 for (
size_t l=0;l<num_local[i];l++) {
677 for (
typename std::set<gno_t>::iterator sitr = adjs2[l].begin();sitr!=adjs2[l].end();sitr++) {
678 adj2_gids[i][j].push_back(*sitr);
680 adj2_offsets[i][j][l+1]=adj2_gids[i][j].size();
687 ent_coords =
new scalar_t*[m_dimension+1];
688 for (
int i=0;i<=m_dimension;i++) {
689 ent_coords[i] = NULL;
692 apf::MeshIterator* itr = m->begin(i);
693 apf::MeshEntity* ent;
694 ent_coords[i] =
new scalar_t[3*num_local[i]];
696 while((ent=m->iterate(itr))) {
699 m->getPoint(ent,0,point);
702 point = apf::getLinearCentroid(m,ent);
704 for (
int k=0;k<3;k++)
705 ent_coords[i][j*3+k] = point[k];
713 weights =
new map_array_t[m_dimension+1];
716 delete [] lid_mapping;
718template <
typename User>
719void APFMeshAdapter<User>::destroy() {
723 for (
int i=0;i<=m_dimension;i++) {
726 delete [] ent_coords[i];
727 delete [] adj_gids[i];
729 delete [] adj2_gids[i];
730 for (
int j=0;j<=m_dimension;j++) {
734 delete [] adj_offsets[i][j];
736 delete [] adj2_offsets[i][j];
740 delete [] adj2_offsets[i];
741 delete [] adj_offsets[i];
742 delete [] gid_mapping[i];
743 apf::destroyGlobalNumbering(gids[i]);
744 apf::destroyNumbering(lids[i]);
746 delete [] ent_coords;
748 delete [] adj_offsets;
751 delete [] adj2_offsets;
753 delete [] gid_mapping;
762template <
typename User>
763void APFMeshAdapter<User>::setWeights(
MeshEntityType etype,
const scalar_t *val,
int stride,
int idx) {
764 int dim = entityZ2toAPF(etype);
765 if (dim>m_dimension||!has(dim)) {
766 throw std::runtime_error(
"Cannot add weights to non existing dimension");
768 ArrayRCP<const scalar_t> weight_rcp(val,0,stride*getLocalNumOf(etype),
false);
769 weights[dim][idx] =std::make_pair(weight_rcp,stride);
773template <
typename User>
774void APFMeshAdapter<User>::getTagWeight(apf::Mesh* m,
776 apf::MeshEntity* ent,
778 int size = m->getTagSize(tag);
779 int type = m->getTagType(tag);
780 if (type==apf::Mesh::DOUBLE) {
781 double* w =
new double[size];
782 m->getDoubleTag(ent,tag,w);
783 for (
int i=0;i<size;i++)
784 ws[i] =
static_cast<scalar_t
>(w[i]);
787 else if (type==apf::Mesh::INT) {
788 int* w =
new int[size];
789 m->getIntTag(ent,tag,w);
790 for (
int i=0;i<size;i++)
791 ws[i] =
static_cast<scalar_t
>(w[i]);
794 else if (type==apf::Mesh::LONG) {
795 long* w =
new long[size];
796 m->getLongTag(ent,tag,w);
797 for (
int i=0;i<size;i++)
798 ws[i] =
static_cast<scalar_t
>(w[i]);
802 throw std::runtime_error(
"Unrecognized tag type");
806template <
typename User>
807void APFMeshAdapter<User>::setWeights(
MeshEntityType etype, apf::Mesh* m,apf::MeshTag* tag,
int* ids) {
808 int dim = entityZ2toAPF(etype);
809 if (dim>m_dimension||!has(dim)) {
810 throw std::runtime_error(
"Cannot add weights to non existing dimension");
812 int n_weights = m->getTagSize(tag);
813 bool delete_ids =
false;
815 ids =
new int[n_weights];
817 for (
int i=0;i<n_weights;i++)
820 scalar_t* ones =
new scalar_t[n_weights];
821 for (
int i=0;i<n_weights;i++)
824 scalar_t* ws =
new scalar_t[num_local[dim]*n_weights];
825 apf::MeshIterator* itr = m->begin(dim);
826 apf::MeshEntity* ent;
828 while ((ent=m->iterate(itr))) {
830 if (m->hasTag(ent,tag)) {
831 w =
new scalar_t[n_weights];
832 getTagWeight(m,tag,ent,w);
837 for (
int i=0;i<n_weights;i++) {
838 ws[i*getLocalNumOf(etype)+j] = w[i];
842 if (m->hasTag(ent,tag))
845 for (
int i=0;i<n_weights;i++) {
846 ArrayRCP<const scalar_t> weight_rcp(ws+i*getLocalNumOf(etype),0,getLocalNumOf(etype),i==0);
847 weights[dim][ids[i]] =std::make_pair(weight_rcp,1);
855template <
typename User>
856void APFMeshAdapter<User>::print(
int me,
int verbosity)
858 if (m_dimension==-1) {
859 std::cout<<
"Cannot print destroyed mesh adapter\n";
863 std::string fn(
" APFMesh ");
864 std::cout << me << fn
865 <<
" dimension = " << m_dimension
869 for (
int i=0;i<=m_dimension;i++) {
872 std::cout<<me<<
" Number of dimension " << i<<
" = " <<num_local[i] <<std::endl;
874 for (
size_t j=0;j<num_local[i];j++) {
875 std::cout<<
" Entity "<<gid_mapping[i][j]<<
"("<<j<<
"):\n";
876 for (
int k=0;k<=m_dimension;k++) {
881 std::cout<<
" First Adjacency of Dimension "<<k<<
":";
882 for (offset_t l=adj_offsets[i][k][j];l<adj_offsets[i][k][j+1];l++)
883 std::cout<<
" "<<adj_gids[i][k][l];
886 std::cout<<
" Second Adjacency through Dimension "<<k<<
":";
887 for (offset_t l=adj2_offsets[i][k][j];l<adj2_offsets[i][k][j+1];l++)
888 std::cout<<
" "<<adj2_gids[i][k][l];
Zoltan2::BasicUserTypes< zscalar_t, zlno_t, zgno_t > user_t
Defines the MeshAdapter interface.
This file defines the StridedData class.
APFMeshAdapter(const Comm< int > &comm, apf::Mesh *m, std::string primary, std::string adjacency, bool needSecondAdj=false)
InputTraits< User >::node_t node_t
InputTraits< User >::offset_t offset_t
InputTraits< User >::part_t part_t
InputTraits< User >::scalar_t scalar_t
InputTraits< User >::lno_t lno_t
void applyPartitioningSolution(const User &in, User *&out, const PartitioningSolution< Adapter > &solution) const
Apply a PartitioningSolution to an input.
InputTraits< User >::gno_t gno_t
MeshAdapter defines the interface for mesh input.
virtual void getWeightsViewOf(MeshEntityType etype, const scalar_t *&weights, int &stride, int idx=0) const
Provide a pointer to one of the number of this process' optional entity weights.
virtual bool availAdjs(MeshEntityType source, MeshEntityType target) const
Returns whether a first adjacency combination is available.
virtual int getNumWeightsPerOf(MeshEntityType etype) const
Return the number of weights per entity.
virtual size_t getLocalNum2ndAdjs(MeshEntityType sourcetarget, MeshEntityType through) const
if avail2ndAdjs(), returns the number of second adjacencies on this process.
virtual void getTopologyViewOf(MeshEntityType etype, enum EntityTopologyType const *&Types) const
Provide a pointer to the entity topology types.
virtual int getDimension() const
Return dimension of the entity coordinates, if any.
virtual size_t getLocalNumOf(MeshEntityType etype) const =0
Returns the global number of mesh entities of MeshEntityType.
virtual size_t getLocalNumAdjs(MeshEntityType source, MeshEntityType target) const
Returns the number of first adjacencies on this process.
virtual void getIDsViewOf(MeshEntityType etype, gno_t const *&Ids) const =0
Provide a pointer to this process' identifiers.
virtual bool areEntityIDsUnique(MeshEntityType etype) const
Provide a pointer to the entity topology types.
virtual void getCoordinatesViewOf(MeshEntityType etype, const scalar_t *&coords, int &stride, int coordDim) const
Provide a pointer to one dimension of entity coordinates.
virtual void getAdjsView(MeshEntityType source, MeshEntityType target, const offset_t *&offsets, const gno_t *&adjacencyIds) const
Sets pointers to this process' mesh first adjacencies.
virtual bool avail2ndAdjs(MeshEntityType sourcetarget, MeshEntityType through) const
Returns whether a second adjacency combination is available. If combination is not available in the M...
enum MeshEntityType getPrimaryEntityType() const
Returns the entity to be partitioned, ordered, colored, etc.
virtual void get2ndAdjsView(MeshEntityType sourcetarget, MeshEntityType through, const offset_t *&offsets, const gno_t *&adjacencyIds) const
if avail2ndAdjs(), set pointers to this process' second adjacencies
map_t::local_ordinal_type lno_t
map_t::global_ordinal_type gno_t
Created by mbenlioglu on Aug 31, 2020.
EntityTopologyType
Enumerate entity topology types for meshes: points,lines,polygons,triangles,quadrilaterals,...
MeshEntityType
Enumerate entity types for meshes: Regions, Faces, Edges, or Vertices.
SparseMatrixAdapter_t::part_t part_t