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> {
142 APFMeshAdapter(
const Comm<int> &comm, apf::Mesh* m,std::string primary,
143 std::string adjacency,
bool needSecondAdj=
false,
int needs=0);
146 void print(
int me,
int verbosity=0);
147 template <
typename Adapter>
149 const PartitioningSolution<Adapter> &solution)
const{
151 apf::Migration* plan =
new apf::Migration(*out);
152 const part_t* new_part_ids = solution.getPartListView();
157 apf::MeshIterator* itr = (*out)->begin(m_dimension);
158 apf::MeshEntity* ent;
160 while ((ent=(*out)->iterate(itr))) {
161 assert(new_part_ids[i]<PCU_Comm_Peers());
162 plan->send(ent,new_part_ids[i]);
171 apf::MeshIterator* itr = (*out)->begin(m_dimension);
172 apf::MeshEntity* ent;
174 while ((ent=(*out)->iterate(itr))) {
175 std::unordered_map<unsigned int,unsigned int> newOwners;
177 unsigned int max_num = 0;
178 int new_part=PCU_Comm_Self();
179 unsigned int num = in->getDownward(ent,dim,adj);
180 for (
unsigned int j=0;j<num;j++) {
182 lno_t lid = apf::getNumber(lids[dim],adj[j],0,0);
183 newOwners[new_part_ids[lid]]++;
184 if (newOwners[new_part_ids[lid]]>max_num) {
185 max_num=newOwners[new_part_ids[lid]];
186 new_part = new_part_ids[lid];
190 if (new_part<0||new_part>=PCU_Comm_Peers()) {
191 std::cout<<new_part<<std::endl;
192 throw std::runtime_error(
"Target part is out of bounds\n");
194 plan->send(ent,new_part);
199 (*out)->migrate(plan);
213 int dim = entityZ2toAPF(etype);
214 return dim==m_dimension;
218 int dim = entityZ2toAPF(etype);
219 if (dim<=m_dimension&&dim>=0)
220 return num_local[dim];
226 int dim = entityZ2toAPF(etype);
227 if (dim<=m_dimension&&dim>=0)
228 Ids = gid_mapping[dim];
235 int dim = entityZ2toAPF(etype);
236 if (dim<=m_dimension&&dim>=0)
237 Types = topologies[dim];
243 int dim = entityZ2toAPF(etype);
244 return static_cast<int>(
weights[dim].size());
248 int &stride,
int idx = 0)
const 250 int dim = entityZ2toAPF(etype);
251 typename map_array_t::iterator itr =
weights[dim].find(idx);
253 ws = &(*(itr->second.first));
254 stride = itr->second.second;
265 int &stride,
int coordDim)
const {
266 if (coordDim>=0 && coordDim<3) {
267 int dim = entityZ2toAPF(etype);
268 if (dim<=m_dimension&&dim>=0) {
269 coords = ent_coords[dim]+coordDim;
284 int dim_source = entityZ2toAPF(source);
285 int dim_target = entityZ2toAPF(target);
286 return dim_source<=m_dimension && dim_source>=0 &&
287 dim_target<=m_dimension && dim_target>=0 &&
288 dim_target!=dim_source&&
289 has(dim_source) && has(dim_target);
294 int dim_source = entityZ2toAPF(source);
295 int dim_target = entityZ2toAPF(target);
297 return adj_gids[dim_source][dim_target].size();
302 const lno_t *&offsets,
const gno_t *& adjacencyIds)
const 304 int dim_source = entityZ2toAPF(source);
305 int dim_target = entityZ2toAPF(target);
307 offsets = adj_offsets[dim_source][dim_target];
308 adjacencyIds = &(adj_gids[dim_source][dim_target][0]);
318 #ifndef USE_MESH_ADAPTER 323 int dim_source = entityZ2toAPF(sourcetarget);
324 int dim_target = entityZ2toAPF(through);
325 if (dim_source==1&&dim_target==0)
327 return dim_source<=m_dimension && dim_source>=0 &&
328 dim_target<=m_dimension && dim_target>=0 &&
329 dim_target!=dim_source &&
330 has(dim_source)&&has(dim_target);
336 int dim_source = entityZ2toAPF(sourcetarget);
337 int dim_target = entityZ2toAPF(through);
339 return adj2_gids[dim_source][dim_target].size();
345 const lno_t *&offsets,
const gno_t *&adjacencyIds)
const 347 int dim_source = entityZ2toAPF(sourcetarget);
348 int dim_target = entityZ2toAPF(through);
350 offsets=adj2_offsets[dim_source][dim_target];
351 adjacencyIds=&(adj2_gids[dim_source][dim_target][0]);
390 bool has(
int dim)
const {
return (entity_needs>>dim)%2;}
393 int entityZ2toAPF(
enum MeshEntityType etype)
const {
return static_cast<int>(etype);}
397 if (ttype==apf::Mesh::VERTEX)
399 else if (ttype==apf::Mesh::EDGE)
403 else if (ttype==apf::Mesh::QUAD)
405 else if (ttype==apf::Mesh::TET)
407 else if (ttype==apf::Mesh::HEX)
414 throw "No such APF topology type";
419 void getTagWeight(apf::Mesh* m, apf::MeshTag* tag,apf::MeshEntity* ent,
scalar_t* ws);
428 apf::Numbering** lids;
429 apf::GlobalNumbering** gids;
433 lno_t*** adj_offsets;
434 std::vector<gno_t>** adj_gids;
435 lno_t*** adj2_offsets;
436 std::vector<gno_t>** adj2_gids;
441 typedef std::unordered_map<int, std::pair<ArrayRCP<const scalar_t>,
int> > map_array_t;
450 template <
typename User>
454 std::string adjacency,
459 m_dimension = m->getDimension();
462 entity_needs = needs;
468 if (primary==
"element") {
474 if (adjacency==
"element") {
480 if (primary==
"region"&&m_dimension<3)
481 throw std::runtime_error(
"primary type and mesh dimension mismatch");
482 if (adjacency==
"region"&&m_dimension<3)
483 throw std::runtime_error(
"adjacency type and mesh dimension mismatch");
484 this->setEntityTypes(primary,adjacency,adjacency);
487 int dim1 = entityZ2toAPF(this->getPrimaryEntityType());
488 int dim2 = entityZ2toAPF(this->getAdjacencyEntityType());
492 entity_needs|=new_needs;
495 lids =
new apf::Numbering*[m_dimension+1];
496 gids =
new apf::GlobalNumbering*[m_dimension+1];
497 gid_mapping =
new gno_t*[m_dimension+1];
498 std::unordered_map<gno_t,lno_t>* lid_mapping =
new std::unordered_map<gno_t,lno_t>[m_dimension+1];
499 num_local =
new size_t[m_dimension+1];
502 for (
int i=0;i<=m_dimension;i++) {
505 topologies[i] = NULL;
506 gid_mapping[i] = NULL;
510 num_local[i] = m->count(i);
511 long global_count = countOwned(m,i);
512 PCU_Add_Longs(&global_count,1);
517 sprintf(lids_name,
"lids%d",i);
519 sprintf(gids_name,
"ids%d",i);
520 apf::FieldShape*
shape = apf::getConstant(i);
521 lids[i] = apf::createNumbering(m,lids_name,
shape,1);
522 apf::Numbering* tmp = apf::numberOwnedDimension(m,gids_name,i);
523 gids[i] = apf::makeGlobal(tmp);
524 apf::synchronize(gids[i]);
525 apf::MeshIterator* itr = m->begin(i);
526 apf::MeshEntity* ent;
528 while ((ent=m->iterate(itr))) {
529 apf::number(lids[i],ent,0,0,num);
530 lid_mapping[i][apf::getNumber(gids[i],
apf::Node(ent,0))]=num;
534 assert(num==num_local[i]);
538 gid_mapping[i] =
new gno_t[num_local[i]];
540 apf::DynamicArray<apf::Node> nodes;
543 while((ent=m->iterate(itr))) {
544 gno_t gid = apf::getNumber(gids[i],
apf::Node(ent,0));
545 gid_mapping[i][ apf::getNumber(lids[i],ent,0,0)] = gid;
554 adj_gids =
new std::vector<gno_t>*[m_dimension+1];
555 adj_offsets =
new lno_t**[m_dimension+1];
557 adj2_gids =
new std::vector<gno_t>*[m_dimension+1];
558 adj2_offsets =
new lno_t**[m_dimension+1];
564 for (
int i=0;i<=m_dimension;i++) {
569 adj2_offsets[i]=NULL;
573 adj_gids[i] =
new std::vector<gno_t>[m_dimension+1];
574 adj_offsets[i] =
new lno_t*[m_dimension+1];
576 adj2_gids[i] =
new std::vector<gno_t>[m_dimension+1];
577 adj2_offsets[i] =
new lno_t*[m_dimension+1];
579 for (
int j=0;j<=m_dimension;j++) {
582 adj_offsets[i][j]=NULL;
584 adj2_offsets[i][j]=NULL;
589 apf::MeshIterator* itr = m->begin(i);
590 apf::MeshEntity* ent;
592 adj_offsets[i][j] =
new lno_t[num_local[i]+1];
593 adj_offsets[i][j][0] =0;
595 adj2_offsets[i][j] =
new lno_t[num_local[i]+1];
596 adj2_offsets[i][j][0] =0;
603 std::unordered_map<gno_t,apf::MeshEntity*> part_boundary_mapping;
605 while ((ent=m->iterate(itr))) {
606 std::set<gno_t> temp_adjs;
609 m->getAdjacent(ent,j,adj);
610 for (
unsigned int l=0;l<adj.getSize();l++) {
611 adj_gids[i][j].push_back(apf::getNumber(gids[j],
apf::Node(adj[l],0)));
615 m->getAdjacent(adj[l],i,adj2);
616 for (
unsigned int o=0;o<adj2.getSize();o++)
617 temp_adjs.insert(apf::getNumber(gids[i],
apf::Node(adj2[o],0)));
618 if (i==m_dimension) {
620 m->getResidence(adj[l],res);
622 part_boundary_mapping[apf::getNumber(gids[j],
apf::Node(adj[l],0))] = adj[l];
623 for (apf::Parts::iterator it=res.begin();it!=res.end();it++) {
625 send_vals[1]=apf::getNumber(gids[i],
apf::Node(ent,0));
626 send_vals[0]=apf::getNumber(gids[j],
apf::Node(adj[l],0));
628 PCU_Comm_Pack(*it,send_vals,2*
sizeof(gno_t));
633 adj_offsets[i][j][k] = adj_gids[i][j].size();
636 if (needSecondAdj && i!=m_dimension) {
638 m->getResidence(ent,res);
639 typename std::set<gno_t>::iterator adj_itr;
640 for (adj_itr=temp_adjs.begin();adj_itr!=temp_adjs.end();adj_itr++) {
641 for (apf::Parts::iterator it=res.begin();it!=res.end();it++) {
643 send_vals[0]=apf::getNumber(gids[i],
apf::Node(ent,0));
644 send_vals[1] = *adj_itr;
645 if (send_vals[0]!=send_vals[1])
646 PCU_Comm_Pack(*it,send_vals,2*
sizeof(gno_t));
655 std::set<gno_t>* adjs2 =
new std::set<gno_t>[num_local[i]];
656 while (PCU_Comm_Receive()) {
658 PCU_Comm_Unpack(adj2,2*
sizeof(gno_t));
660 if (i==m_dimension) {
661 apf::MeshEntity* through = part_boundary_mapping[adj2[0]];
663 m->getAdjacent(through,i,adj);
664 for (
unsigned int l=0;l<adj.getSize();l++) {
665 if (apf::getNumber(gids[i],
apf::Node(adj[l],0))!=adj2[1])
666 adjs2[apf::getNumber(lids[i],adj[l],0,0)].insert(adj2[1]);
670 lno_t index = lid_mapping[i][adj2[0]];
671 adjs2[index].insert(adj2[1]);
675 for (
size_t l=0;l<num_local[i];l++) {
676 for (
typename std::set<gno_t>::iterator sitr = adjs2[l].begin();sitr!=adjs2[l].end();sitr++) {
677 adj2_gids[i][j].push_back(*sitr);
679 adj2_offsets[i][j][l+1]=adj2_gids[i][j].size();
686 ent_coords =
new scalar_t*[m_dimension+1];
687 for (
int i=0;i<=m_dimension;i++) {
688 ent_coords[i] = NULL;
691 apf::MeshIterator* itr = m->begin(i);
692 apf::MeshEntity* ent;
693 ent_coords[i] =
new scalar_t[3*num_local[i]];
695 while((ent=m->iterate(itr))) {
698 m->getPoint(ent,0,point);
701 point = apf::getLinearCentroid(m,ent);
703 for (
int k=0;k<3;k++)
704 ent_coords[i][j*3+k] = point[k];
712 weights =
new map_array_t[m_dimension+1];
715 delete [] lid_mapping;
717 template <
typename User>
718 void APFMeshAdapter<User>::destroy() {
722 for (
int i=0;i<=m_dimension;i++) {
725 delete [] ent_coords[i];
726 delete [] adj_gids[i];
728 delete [] adj2_gids[i];
729 for (
int j=0;j<=m_dimension;j++) {
733 delete [] adj_offsets[i][j];
735 delete [] adj2_offsets[i][j];
739 delete [] adj2_offsets[i];
740 delete [] adj_offsets[i];
741 delete [] gid_mapping[i];
742 apf::destroyGlobalNumbering(gids[i]);
743 apf::destroyNumbering(lids[i]);
745 delete [] ent_coords;
747 delete [] adj_offsets;
750 delete [] adj2_offsets;
752 delete [] gid_mapping;
761 template <
typename User>
762 void APFMeshAdapter<User>::setWeights(
MeshEntityType etype,
const scalar_t *val,
int stride,
int idx) {
763 int dim = entityZ2toAPF(etype);
764 if (dim>m_dimension||!has(dim)) {
765 throw std::runtime_error(
"Cannot add weights to non existing dimension");
767 ArrayRCP<const scalar_t> weight_rcp(val,0,stride*getLocalNumOf(etype),
false);
768 weights[dim][idx] =std::make_pair(weight_rcp,stride);
772 template <
typename User>
773 void APFMeshAdapter<User>::getTagWeight(apf::Mesh* m,
775 apf::MeshEntity* ent,
777 int size = m->getTagSize(tag);
778 int type = m->getTagType(tag);
779 if (type==apf::Mesh::DOUBLE) {
780 double* w =
new double[size];
781 m->getDoubleTag(ent,tag,w);
782 for (
int i=0;i<size;i++)
783 ws[i] = static_cast<scalar_t>(w[i]);
786 else if (type==apf::Mesh::INT) {
787 int* w =
new int[size];
788 m->getIntTag(ent,tag,w);
789 for (
int i=0;i<size;i++)
790 ws[i] = static_cast<scalar_t>(w[i]);
793 else if (type==apf::Mesh::LONG) {
794 long* w =
new long[size];
795 m->getLongTag(ent,tag,w);
796 for (
int i=0;i<size;i++)
797 ws[i] = static_cast<scalar_t>(w[i]);
801 throw std::runtime_error(
"Unrecognized tag type");
805 template <
typename User>
806 void APFMeshAdapter<User>::setWeights(
MeshEntityType etype, apf::Mesh* m,apf::MeshTag* tag,
int* ids) {
807 int dim = entityZ2toAPF(etype);
808 if (dim>m_dimension||!has(dim)) {
809 throw std::runtime_error(
"Cannot add weights to non existing dimension");
811 int n_weights = m->getTagSize(tag);
812 bool delete_ids =
false;
814 ids =
new int[n_weights];
816 for (
int i=0;i<n_weights;i++)
819 scalar_t* ones =
new scalar_t[n_weights];
820 for (
int i=0;i<n_weights;i++)
823 scalar_t* ws =
new scalar_t[num_local[dim]*n_weights];
824 apf::MeshIterator* itr = m->begin(dim);
825 apf::MeshEntity* ent;
827 while ((ent=m->iterate(itr))) {
829 if (m->hasTag(ent,tag)) {
830 w =
new scalar_t[n_weights];
831 getTagWeight(m,tag,ent,w);
836 for (
int i=0;i<n_weights;i++) {
837 ws[i*getLocalNumOf(etype)+j] = w[i];
841 if (m->hasTag(ent,tag))
844 for (
int i=0;i<n_weights;i++) {
845 ArrayRCP<const scalar_t> weight_rcp(ws+i*getLocalNumOf(etype),0,getLocalNumOf(etype),i==0);
846 weights[dim][ids[i]] =std::make_pair(weight_rcp,1);
854 template <
typename User>
855 void APFMeshAdapter<User>::print(
int me,
int verbosity)
857 if (m_dimension==-1) {
858 std::cout<<
"Cannot print destroyed mesh adapter\n";
862 std::string fn(
" APFMesh ");
863 std::cout << me << fn
864 <<
" dimension = " << m_dimension
868 for (
int i=0;i<=m_dimension;i++) {
871 std::cout<<me<<
" Number of dimension " << i<<
" = " <<num_local[i] <<std::endl;
873 for (
size_t j=0;j<num_local[i];j++) {
874 std::cout<<
" Entity "<<gid_mapping[i][j]<<
"("<<j<<
"):\n";
875 for (
int k=0;k<=m_dimension;k++) {
880 std::cout<<
" First Adjacency of Dimension "<<k<<
":";
881 for (lno_t l=adj_offsets[i][k][j];l<adj_offsets[i][k][j+1];l++)
882 std::cout<<
" "<<adj_gids[i][k][l];
885 std::cout<<
" Second Adjacency through Dimension "<<k<<
":";
886 for (lno_t l=adj2_offsets[i][k][j];l<adj2_offsets[i][k][j+1];l++)
887 std::cout<<
" "<<adj2_gids[i][k][l];
900 #endif //HAVE_ZOLTAN2_PARMA InputTraits< User >::scalar_t scalar_t
virtual int getNumWeightsPerOf(MeshEntityType etype) const
Return the number of weights per entity.
APFMeshAdapter(const Comm< int > &comm, apf::Mesh *m, std::string primary, std::string adjacency, bool needSecondAdj=false)
virtual size_t getLocalNum2ndAdjs(MeshEntityType sourcetarget, MeshEntityType through) const
if avail2ndAdjs(), returns the number of second adjacencies on this process.
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 bool availAdjs(MeshEntityType source, MeshEntityType target) const
Returns whether a first adjacency combination is available.
virtual void getTopologyViewOf(MeshEntityType etype, enum EntityTopologyType const *&Types) const
Provide a pointer to the entity topology types.
InputTraits< User >::lno_t lno_t
virtual size_t getLocalNumAdjs(MeshEntityType source, MeshEntityType target) const
Returns the number of first adjacencies on this process.
void applyPartitioningSolution(const User &in, User *&out, const PartitioningSolution< Adapter > &solution) const
Apply a PartitioningSolution to an input.
virtual void getIDsViewOf(MeshEntityType etype, gno_t const *&Ids) const =0
Provide a pointer to this process' identifiers.
virtual void getAdjsView(MeshEntityType source, MeshEntityType target, const lno_t *&offsets, const gno_t *&adjacencyIds) const
Sets pointers to this process' mesh first adjacencies.
InputTraits< User >::part_t part_t
EntityTopologyType
Enumerate entity topology types for meshes: points,lines,polygons,triangles,quadrilaterals, polyhedrons, tetrahedrons, hexhedrons, prisms, or pyramids.
virtual void get2ndAdjsView(MeshEntityType sourcetarget, MeshEntityType through, const lno_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.
enum MeshEntityType getPrimaryEntityType() const
Returns the entity to be partitioned, ordered, colored, etc.
MeshEntityType
Enumerate entity types for meshes: Regions, Faces, Edges, or Vertices.
virtual void getCoordinatesViewOf(MeshEntityType etype, const scalar_t *&coords, int &stride, int coordDim) const
Provide a pointer to one dimension of entity coordinates.
virtual bool areEntityIDsUnique(MeshEntityType etype) const
Provide a pointer to the entity topology types.
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 avail2ndAdjs(MeshEntityType sourcetarget, MeshEntityType through) const
Returns whether a second adjacency combination is available. If combination is not available in the M...
This file defines the StridedData class.