50 #ifndef _ZOLTAN2_APFMESHADAPTER_HPP_
51 #define _ZOLTAN2_APFMESHADAPTER_HPP_
56 #include <unordered_map>
61 #ifndef HAVE_ZOLTAN2_PARMA
67 template <
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>
118 template <
typename User>
119 class APFMeshAdapter:
public MeshAdapter<User> {
123 typedef typename InputTraits<User>::scalar_t
scalar_t;
124 typedef typename InputTraits<User>::lno_t
lno_t;
125 typedef typename InputTraits<User>::gno_t
gno_t;
126 typedef typename InputTraits<User>::offset_t
offset_t;
128 typedef typename InputTraits<User>::node_t node_t;
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++) {
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)
404 else if (ttype==apf::Mesh::QUAD)
406 else if (ttype==apf::Mesh::TET)
408 else if (ttype==apf::Mesh::HEX)
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;
451 template <
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;
718 template <
typename User>
719 void 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;
762 template <
typename User>
763 void 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);
773 template <
typename User>
774 void 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");
806 template <
typename User>
807 void 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);
855 template <
typename User>
856 void 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];
901 #endif //HAVE_ZOLTAN2_PARMA
virtual void getAdjsView(MeshEntityType source, MeshEntityType target, const offset_t *&offsets, const gno_t *&adjacencyIds) const
Sets pointers to this process' mesh first adjacencies.
InputTraits< User >::scalar_t scalar_t
APFMeshAdapter(const Comm< int > &comm, apf::Mesh *m, std::string primary, std::string adjacency, bool needSecondAdj=false)
virtual bool areEntityIDsUnique(MeshEntityType etype) const
Provide a pointer to the entity topology types.
virtual bool availAdjs(MeshEntityType source, MeshEntityType target) const
Returns whether a first adjacency combination is available.
InputTraits< User >::gno_t gno_t
Defines the MeshAdapter interface.
MeshAdapter defines the interface for mesh input.
virtual size_t getLocalNumOf(MeshEntityType etype) const =0
Returns the global number of mesh entities of MeshEntityType.
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.
InputTraits< User >::lno_t lno_t
virtual bool avail2ndAdjs(MeshEntityType sourcetarget, MeshEntityType through) const
Returns whether a second adjacency combination is available. If combination is not available in the M...
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
virtual int getDimension() const
Return dimension of the entity coordinates, if any.
SparseMatrixAdapter_t::part_t part_t
virtual size_t getLocalNum2ndAdjs(MeshEntityType sourcetarget, MeshEntityType through) const
if avail2ndAdjs(), returns the number of second adjacencies on this process.
virtual void getIDsViewOf(MeshEntityType etype, gno_t const *&Ids) const =0
Provide a pointer to this process' identifiers.
virtual void getCoordinatesViewOf(MeshEntityType etype, const scalar_t *&coords, int &stride, int coordDim) const
Provide a pointer to one dimension of entity coordinates.
virtual size_t getLocalNumAdjs(MeshEntityType source, MeshEntityType target) const
Returns the number of first adjacencies on this process.
InputTraits< User >::part_t part_t
EntityTopologyType
Enumerate entity topology types for meshes: points,lines,polygons,triangles,quadrilaterals, polyhedrons, tetrahedrons, hexhedrons, prisms, or pyramids.
void applyPartitioningSolution(const User &in, User *&out, const PartitioningSolution< Adapter > &solution) const
Apply a PartitioningSolution to an input.
MeshEntityType
Enumerate entity types for meshes: Regions, Faces, Edges, or Vertices.
virtual int getNumWeightsPerOf(MeshEntityType etype) const
Return the number of weights per entity.
virtual void getTopologyViewOf(MeshEntityType etype, enum EntityTopologyType const *&Types) const
Provide a pointer to the entity topology types.
enum MeshEntityType getPrimaryEntityType() const
Returns the entity to be partitioned, ordered, colored, etc.
InputTraits< User >::offset_t offset_t
This file defines the StridedData class.
Zoltan2::BasicUserTypes< zscalar_t, zlno_t, zgno_t > user_t