Actual source code: Completion.hh
petsc-3.4.2 2013-07-02

      1: #ifndef included_ALE_Completion_hh
      2: #define included_ALE_Completion_hh

      4: #ifndef  included_ALE_Sections_hh
      5: #include <sieve/Sections.hh>
      6: #endif

      8: #ifndef  included_ALE_SectionCompletion_hh
      9: #include <sieve/SectionCompletion.hh>
     10: #endif

     12: #ifndef  included_ALE_ParallelMapping_hh
     13: #include <sieve/ParallelMapping.hh>
     14: #endif

     16: #include <iostream>
     17: #include <fstream>

     19: namespace ALE {
     20:   class Completion {
     21:   public:
     22:     template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
     23:     static void completeSection(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
     24:       typedef ALE::Section<ALE::Pair<int, typename SendOverlap::source_type>, typename SendSection::value_type> OverlapSection;
     25:       //typedef ALE::Section<typename SendSection::point_type, typename SendSection::value_type> OverlapSection;
     26:       Obj<OverlapSection> overlapSection = new OverlapSection(sendSection->comm(), sendSection->debug());

     28:       if (sendSection->debug()) {sendSection->view("Send Section");}
     29:       ALE::Pullback::SimpleCopy::copy(sendOverlap, recvOverlap, sendSection, overlapSection);
     30:       if (overlapSection->debug()) {overlapSection->view("Overlap Section");}
     31:       ALE::Pullback::InsertionBinaryFusion::fuse(overlapSection, recvOverlap, recvSection);
     32:       if (recvSection->debug()) {recvSection->view("Receieve Section");}
     33:     }
     34:     template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
     35:     static void completeSectionAdd(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
     36:       typedef ALE::Section<ALE::Pair<int, typename SendOverlap::source_type>, typename SendSection::value_type> OverlapSection;
     37:       Obj<OverlapSection> overlapSection = new OverlapSection(sendSection->comm(), sendSection->debug());

     39:       if (sendSection->debug()) {sendSection->view("Send Section");}
     40:       ALE::Pullback::SimpleCopy::copy(sendOverlap, recvOverlap, sendSection, overlapSection);
     41:       if (overlapSection->debug()) {overlapSection->view("Overlap Section");}
     42:       ALE::Pullback::AdditiveBinaryFusion::fuse(overlapSection, recvOverlap, recvSection);
     43:       if (recvSection->debug()) {recvSection->view("Receieve Section");}
     44:     }
     45:     template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection, typename OverlapSection>
     46:     static void completeSectionAdd(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection, const Obj<OverlapSection>& overlapSection) {
     47:       if (sendSection->debug()) {sendSection->view("Send Section");}
     48:       ALE::Pullback::SimpleCopy::copy(sendOverlap, recvOverlap, sendSection, overlapSection);
     49:       if (overlapSection->debug()) {overlapSection->view("Overlap Section");}
     50:       ALE::Pullback::AdditiveBinaryFusion::fuse(overlapSection, recvOverlap, recvSection);
     51:       if (recvSection->debug()) {recvSection->view("Receieve Section");}
     52:     }
     53:   };
     54:   namespace New {
     55:     template<typename Bundle_, typename Value_, typename Alloc_ = malloc_allocator<typename Bundle_::point_type> >
     56:     class Completion {
     57:     public:
     58:       typedef int                                                                         point_type;
     59:       typedef Value_                                                                      value_type;
     60:       typedef Bundle_                                                                     bundle_type;
     61:       typedef Alloc_                                                                      alloc_type;
     62:       typedef typename alloc_type::template rebind<int>::other                            int_alloc_type;
     63:       typedef typename alloc_type::template rebind<value_type>::other                     value_alloc_type;
     64:       typedef typename bundle_type::sieve_type                                            sieve_type;
     65:       typedef typename ALE::DiscreteSieve<point_type, alloc_type>                         dsieve_type;
     66:       typedef typename ALE::Topology<int, dsieve_type, alloc_type>                        topology_type;
     67:       typedef typename bundle_type::rank_type                                             rank_type;
     68:       typedef typename ALE::Sifter<point_type, rank_type, point_type>                     send_overlap_type;
     69:       typedef typename ALE::Sifter<rank_type, point_type, point_type>                     recv_overlap_type;
     70:       typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, int, int_alloc_type> > send_sizer_type;
     71:       typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, int, int_alloc_type> > recv_sizer_type;
     72:       typedef typename ALE::New::ConeSizeSection<bundle_type, sieve_type>                 cone_size_section;
     73:       typedef typename ALE::New::ConeSection<sieve_type>                                  cone_section;
     74:       typedef typename ALE::New::SectionCompletion<bundle_type, value_type, alloc_type>   completion;
     75:     public:
     76:       template<typename PartitionType>
     77:       static void scatterSieve(const Obj<bundle_type>& bundle, const Obj<sieve_type>& sieve, const int dim, const Obj<sieve_type>& sieveNew, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const int height, const int numCells, PartitionType assignment[]) {
     78:         typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > send_section_type;
     79:         typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > recv_section_type;
     80:         int rank  = sieve->commRank();
     81:         int debug = sieve->debug();

     83:         // Create local sieve
     84:         const Obj<typename bundle_type::label_sequence>& cells = bundle->heightStratum(height);
     85:         int e = 0;

     87:         if (sieve->debug()) {
     88:           int e2 = 0;
     89:           for(typename bundle_type::label_sequence::iterator e_iter = cells->begin(); e_iter != cells->end(); ++e_iter) {
     90:             std::cout << "assignment["<<*e_iter<<"]" << assignment[e2++] << std::endl;
     91:           }
     92:         }
     93:         PetscBool  flg;
     94:         PetscOptionsHasName(NULL, "-output_partition", &flg);
     95:         if (flg) {
     96:           ostringstream fname;
     97:           fname << "part." << sieve->commSize() << ".dat";
     98:           std::ofstream f(fname.str().c_str());
     99:           int e2 = 0;
    100:           f << sieve->commSize() << std::endl;
    101:           for(typename bundle_type::label_sequence::iterator e_iter = cells->begin(); e_iter != cells->end(); ++e_iter) {
    102:             f << assignment[e2++] << std::endl;
    103:           }
    104:         }
    105:         for(typename bundle_type::label_sequence::iterator e_iter = cells->begin(); e_iter != cells->end(); ++e_iter) {
    106:           if (assignment[e] == rank) {
    107:             Obj<typename sieve_type::coneSet> current = new typename sieve_type::coneSet();
    108:             Obj<typename sieve_type::coneSet> next    = new typename sieve_type::coneSet();
    109:             Obj<typename sieve_type::coneSet> tmp;

    111:             current->insert(*e_iter);
    112:             while(current->size()) {
    113:               for(typename sieve_type::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
    114:                 const Obj<typename sieve_type::traits::coneSequence>& cone = sieve->cone(*p_iter);

    116:                 for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
    117:                   sieveNew->addArrow(*c_iter, *p_iter, c_iter.color());
    118:                   next->insert(*c_iter);
    119:                 }
    120:               }
    121:               tmp = current; current = next; next = tmp;
    122:               next->clear();
    123:             }
    124:             if (height) {
    125:               current->insert(*e_iter);
    126:               while(current->size()) {
    127:                 for(typename sieve_type::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
    128:                   const Obj<typename sieve_type::traits::supportSequence>& support = sieve->support(*p_iter);

    130:                   for(typename sieve_type::traits::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter) {
    131:                     sieveNew->addArrow(*p_iter, *s_iter, s_iter.color());
    132:                     next->insert(*s_iter);
    133:                   }
    134:                 }
    135:                 tmp = current; current = next; next = tmp;
    136:                 next->clear();
    137:               }
    138:             }
    139:           }
    140:           e++;
    141:         }
    142:         // Complete sizer section
    143:         typedef typename ALE::New::PartitionSizeSection<bundle_type, PartitionType> partition_size_section;
    144:         typedef typename ALE::New::PartitionSection<bundle_type, PartitionType>     partition_section;
    145:         Obj<topology_type>          secTopology          = completion::createSendTopology(sendOverlap);
    146:         Obj<partition_size_section> partitionSizeSection = new partition_size_section(bundle, height, numCells, assignment);
    147:         Obj<partition_section>      partitionSection     = new partition_section(bundle, height, numCells, assignment);
    148:         Obj<send_section_type>      sendSection          = new send_section_type(sieve->comm(), sieve->debug());
    149:         Obj<recv_section_type>      recvSection          = new recv_section_type(sieve->comm(), sendSection->getTag(), sieve->debug());

    151:         completion::completeSection(sendOverlap, recvOverlap, partitionSizeSection, partitionSection, sendSection, recvSection);
    152:         // Unpack the section into the overlap
    153:         sendOverlap->clear();
    154:         recvOverlap->clear();
    155:         const typename send_section_type::sheaf_type& sendPatches = sendSection->getPatches();

    157:         for(typename send_section_type::sheaf_type::const_iterator p_iter = sendPatches.begin(); p_iter != sendPatches.end(); ++p_iter) {
    158:           const typename send_section_type::patch_type               rank    = p_iter->first;
    159:           const Obj<typename send_section_type::section_type>&       section = p_iter->second;
    160:           const typename send_section_type::section_type::chart_type chart   = section->getChart();

    162:           for(typename send_section_type::section_type::chart_type::iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
    163:             const typename send_section_type::value_type *points = section->restrictPoint(*c_iter);
    164:             int                                           size   = section->getFiberDimension(*c_iter);

    166:             for(int p = 0; p < size; p++) {
    167:               sendOverlap->addArrow(points[p], rank, points[p]);
    168:             }
    169:           }
    170:         }
    171:         const typename recv_section_type::sheaf_type& recvPatches = recvSection->getPatches();

    173:         for(typename recv_section_type::sheaf_type::const_iterator p_iter = recvPatches.begin(); p_iter != recvPatches.end(); ++p_iter) {
    174:           const typename send_section_type::patch_type               rank    = p_iter->first;
    175:           const Obj<typename send_section_type::section_type>&       section = p_iter->second;
    176:           const typename send_section_type::section_type::chart_type chart   = section->getChart();

    178:           for(typename recv_section_type::section_type::chart_type::iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
    179:             const typename recv_section_type::value_type *points = section->restrictPoint(*c_iter);
    180:             int                                           size   = section->getFiberDimension(*c_iter);

    182:             for(int p = 0; p < size; p++) {
    183:               recvOverlap->addArrow(rank, points[p], points[p]);
    184:             }
    185:           }
    186:         }
    187:         if (debug) {
    188:           sendOverlap->view(std::cout, "Send overlap for points");
    189:           recvOverlap->view(std::cout, "Receive overlap for points");
    190:         }
    191:         // Receive the point section
    192:         ALE::New::Completion<bundle_type, value_type>::scatterCones(sieve, sieveNew, sendOverlap, recvOverlap, bundle, height);
    193:         if (height) {
    194:           ALE::New::Completion<bundle_type, value_type>::scatterSupports(sieve, sieveNew, sendOverlap, recvOverlap, bundle, bundle->depth()-height);
    195:         }
    196:       }
    197:       template<typename SifterType>
    198:       static void scatterCones(const Obj<SifterType>& sifter, const Obj<SifterType>& sifterNew, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const Obj<bundle_type>& bundle = NULL, const int minimumHeight = 0) {
    199:         typedef typename ALE::New::ConeSizeSection<bundle_type, SifterType> cone_size_section;
    200:         typedef typename ALE::New::ConeSection<SifterType>                  cone_section;
    201:         typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > send_section_type;
    202:         typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > recv_section_type;
    203:         Obj<topology_type>     secTopology     = completion::createSendTopology(sendOverlap);
    204:         Obj<cone_size_section> coneSizeSection = new cone_size_section(bundle, sifter, minimumHeight);
    205:         Obj<cone_section>      coneSection     = new cone_section(sifter);
    206:         Obj<send_section_type> sendSection     = new send_section_type(sifter->comm(), sifter->debug());
    207:         Obj<recv_section_type> recvSection     = new recv_section_type(sifter->comm(), sendSection->getTag(), sifter->debug());

    209:         completion::completeSection(sendOverlap, recvOverlap, coneSizeSection, coneSection, sendSection, recvSection);
    210:         // Unpack the section into the sieve
    211:         const typename recv_section_type::sheaf_type& patches = recvSection->getPatches();

    213:         for(typename recv_section_type::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
    214:           const Obj<typename recv_section_type::section_type>&        section = p_iter->second;
    215:           const typename recv_section_type::section_type::chart_type& chart   = section->getChart();

    217:           for(typename recv_section_type::section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
    218:             const typename recv_section_type::value_type *points = section->restrictPoint(*c_iter);
    219:             int size = section->getFiberDimension(*c_iter);
    220:             int c    = 0;

    222:             for(int p = 0; p < size; p++) {
    223:               sifterNew->addArrow(points[p], *c_iter, c++);
    224:             }
    225:           }
    226:         }
    227:       }
    228:       template<typename SifterType, typename Renumbering>
    229:       static void scatterCones(const Obj<SifterType>& sifter, const Obj<SifterType>& sifterNew, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, Renumbering& renumbering, const Obj<bundle_type>& bundle = NULL) {
    230:         PETSc::Log::Event("ScatterCones").begin();
    231:         typedef typename ALE::New::ConeSizeSection<bundle_type, SifterType> cone_size_section;
    232:         typedef typename ALE::New::ConeSection<SifterType>                  cone_section;
    233:         typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > send_section_type;
    234:         typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > recv_section_type;
    235:         Obj<topology_type>     secTopology     = completion::createSendTopology(sendOverlap);
    236:         Obj<cone_size_section> coneSizeSection = new cone_size_section(bundle, sifter);
    237:         Obj<cone_section>      coneSection     = new cone_section(sifter);
    238:         Obj<send_section_type> sendSection     = new send_section_type(sifter->comm(), sifter->debug());
    239:         Obj<recv_section_type> recvSection     = new recv_section_type(sifter->comm(), sendSection->getTag(), sifter->debug());

    241:         PETSc::Log::Event("ScatterConesComplete").begin();
    242:         completion::completeSection(sendOverlap, recvOverlap, coneSizeSection, coneSection, sendSection, recvSection);
    243:         PETSc::Log::Event("ScatterConesComplete").end();
    244:         PETSc::Log::Event("ScatterConesUpdate").begin();
    245:         // Unpack the section into the sieve
    246:         const typename recv_section_type::sheaf_type& patches = recvSection->getPatches();

    248:         for(typename recv_section_type::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
    249:           const Obj<typename recv_section_type::section_type>&        section = p_iter->second;
    250:           const typename recv_section_type::section_type::chart_type& chart   = section->getChart();

    252:           for(typename recv_section_type::section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
    253:             const typename recv_section_type::value_type *points = section->restrictPoint(*c_iter);
    254:             int size = section->getFiberDimension(*c_iter);
    255:             int c    = 0;

    257:             for(int p = 0; p < size; p++) {
    258:               sifterNew->addArrow(points[p], renumbering[*c_iter], c++);
    259:             }
    260:           }
    261:         }
    262:         PETSc::Log::Event("ScatterConesUpdate").end();
    263:         PETSc::Log::Event("ScatterCones").end();
    264:       }
    265:       template<typename SifterType>
    266:       static void scatterSupports(const Obj<SifterType>& sifter, const Obj<SifterType>& sifterNew, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const Obj<bundle_type>& bundle = NULL, const int minimumDepth = 0) {
    267:         typedef typename ALE::New::SupportSizeSection<bundle_type, SifterType> support_size_section;
    268:         typedef typename ALE::New::SupportSection<SifterType>                  support_section;
    269:         typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > send_section_type;
    270:         typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > recv_section_type;
    271:         Obj<topology_type>        secTopology        = completion::createSendTopology(sendOverlap);
    272:         Obj<support_size_section> supportSizeSection = new support_size_section(bundle, sifter, minimumDepth);
    273:         Obj<support_section>      supportSection     = new support_section(sifter);
    274:         Obj<send_section_type>    sendSection        = new send_section_type(sifter->comm(), sifter->debug());
    275:         Obj<recv_section_type>    recvSection        = new recv_section_type(sifter->comm(), sendSection->getTag(), sifter->debug());

    277:         completion::completeSection(sendOverlap, recvOverlap, supportSizeSection, supportSection, sendSection, recvSection);
    278:         // Unpack the section into the sieve
    279:         const typename recv_section_type::sheaf_type& recvPatches = recvSection->getPatches();

    281:         for(typename recv_section_type::sheaf_type::const_iterator p_iter = recvPatches.begin(); p_iter != recvPatches.end(); ++p_iter) {
    282:           const Obj<typename send_section_type::section_type>&       section = p_iter->second;
    283:           const typename send_section_type::section_type::chart_type chart   = section->getChart();

    285:           for(typename recv_section_type::section_type::chart_type::iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
    286:             const typename recv_section_type::value_type *points = section->restrictPoint(*c_iter);
    287:             int                                           size   = section->getFiberDimension(*c_iter);
    288:             int                                           c      = 0;

    290:             for(int p = 0; p < size; p++) {
    291:               sifterNew->addArrow(*c_iter, points[p], c++);
    292:             }
    293:           }
    294:         }
    295:       }
    296:     };

    298:     template<typename Value_>
    299:     class ParallelFactory {
    300:     public:
    301:       typedef Value_ value_type;
    302:     protected:
    303:       int          _debug;
    304:       MPI_Datatype _mpiType;
    305:     protected:
    306:       MPI_Datatype constructMPIType() {
    307:         if (sizeof(value_type) == 1) {
    308:           return MPI_BYTE;
    309:         } else if (sizeof(value_type) == 2) {
    310:           return MPI_SHORT;
    311:         } else if (sizeof(value_type) == 4) {
    312:           return MPI_INT;
    313:         } else if (sizeof(value_type) == 8) {
    314:           return MPI_DOUBLE;
    315:         } else if (sizeof(value_type) == 28) {
    316:           int          blen[2];
    317:           MPI_Aint     indices[2];
    318:           MPI_Datatype oldtypes[2], newtype;
    319:           blen[0] = 1; indices[0] = 0;           oldtypes[0] = MPI_INT;
    320:           blen[1] = 3; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
    321:           MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
    322:           MPI_Type_commit(&newtype);
    323:           return newtype;
    324:         } else if (sizeof(value_type) == 32) {
    325:           int          blen[2];
    326:           MPI_Aint     indices[2];
    327:           MPI_Datatype oldtypes[2], newtype;
    328:           blen[0] = 1; indices[0] = 0;           oldtypes[0] = MPI_DOUBLE;
    329:           blen[1] = 3; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
    330:           MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
    331:           MPI_Type_commit(&newtype);
    332:           return newtype;
    333:         }
    334:         throw ALE::Exception("Cannot determine MPI type for value type");
    335:       };
    336:       ParallelFactory(const int debug) : _debug(debug) {
    337:         this->_mpiType = this->constructMPIType();
    338:       };
    339:     public:
    340:       ~ParallelFactory() {};
    341:     public:
    342:       static const Obj<ParallelFactory>& singleton(const int debug, bool cleanup = false) {
    343:         static Obj<ParallelFactory> *_singleton = NULL;

    345:         if (cleanup) {
    346:           if (debug) {std::cout << "Destroying ParallelFactory" << std::endl;}
    347:           if (_singleton) {delete _singleton;}
    348:           _singleton = NULL;
    349:         } else if (_singleton == NULL) {
    350:           if (debug) {std::cout << "Creating new ParallelFactory" << std::endl;}
    351:           _singleton  = new Obj<ParallelFactory>();
    352:           *_singleton = new ParallelFactory(debug);
    353:         }
    354:         return *_singleton;
    355:       };
    356:     public: // Accessors
    357:       int debug() const {return this->_debug;};
    358:       MPI_Datatype getMPIType() const {return this->_mpiType;};
    359:     };
    360:   }
    361: }
    362: #endif