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

      1: #ifndef included_ALE_SectionCompletion_hh
      2: #define included_ALE_SectionCompletion_hh

      4: #ifndef  included_ALE_Topology_hh
      5: #include <sieve/Topology.hh>
      6: #endif

      8: #ifndef  included_ALE_Field_hh
      9: #include <sieve/Field.hh>
     10: #endif

     12: namespace ALE {
     13:   namespace New {
     14:     template<typename Topology_, typename Value_, typename Alloc_ = typename Topology_::alloc_type>
     15:     class SectionCompletion {
     16:     public:
     17:       typedef int                                                                              point_type;
     18:       typedef Topology_                                                                        mesh_topology_type;
     19:       typedef Value_                                                                           value_type;
     20:       typedef Alloc_                                                                           alloc_type;
     21:       typedef typename alloc_type::template rebind<point_type>::other                          point_alloc_type;
     22:       typedef typename mesh_topology_type::sieve_type                                          sieve_type;
     23:       typedef typename ALE::DiscreteSieve<point_type, point_alloc_type>                        dsieve_type;
     24:       typedef typename ALE::Topology<int, dsieve_type, alloc_type>                             topology_type;
     25:       // TODO: Fix this typedef typename ALE::Partitioner<>::part_type                         rank_type;
     26:       typedef short int                                                                        rank_type;
     27:       typedef typename ALE::New::SectionCompletion<mesh_topology_type, int, alloc_type>        int_completion;
     28:       typedef typename ALE::New::SectionCompletion<mesh_topology_type, value_type, alloc_type> completion;
     29:     public:
     30:       // Creates a DiscreteTopology with the overlap information
     31:       template<typename SendOverlap>
     32:       static Obj<topology_type> createSendTopology(const Obj<SendOverlap>& sendOverlap) {
     33:         const Obj<typename SendOverlap::baseSequence> ranks = sendOverlap->base();
     34:         Obj<topology_type> topology = new topology_type(sendOverlap->comm(), sendOverlap->debug());

     36:         for(typename SendOverlap::baseSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
     37:           Obj<dsieve_type> sendSieve = new dsieve_type(sendOverlap->cone(*r_iter));
     38:           topology->setPatch(*r_iter, sendSieve);
     39:         }
     40:         return topology;
     41:       }
     42:       template<typename RecvOverlap>
     43:       static Obj<topology_type> createRecvTopology(const Obj<RecvOverlap>& recvOverlap) {
     44:         const Obj<typename RecvOverlap::capSequence> ranks = recvOverlap->cap();
     45:         Obj<topology_type> topology = new topology_type(recvOverlap->comm(), recvOverlap->debug());

     47:         for(typename RecvOverlap::capSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
     48:           Obj<dsieve_type> recvSieve = new dsieve_type();
     49:           const Obj<typename RecvOverlap::supportSequence>& points  = recvOverlap->support(*r_iter);

     51:           for(typename RecvOverlap::supportSequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
     52:             recvSieve->addPoint(p_iter.color());
     53:           }
     54:           topology->setPatch(*r_iter, recvSieve);
     55:         }
     56:         return topology;
     57:       }
     58:       template<typename Sizer, typename Section>
     59:       static void setupSend(const Obj<topology_type>& sendChart, const Obj<Sizer>& sendSizer, const Obj<Section>& sendSection) {
     60:         // Here we should just use the overlap as the topology (once it is a new-style sieve)
     61:         sendSection->clear();
     62:         sendSection->setTopology(sendChart);
     63:         sendSection->construct(sendSizer);
     64:         sendSection->allocate();
     65:         if (sendSection->debug() > 10) {sendSection->view("Send section after setup", MPI_COMM_SELF);}
     66:         sendSection->constructCommunication(Section::SEND);
     67:       }
     68:       template<typename Filler, typename Section>
     69:       static void fillSend(const Filler& sendFiller, const Obj<Section>& sendSection) {
     70:         const typename Section::sheaf_type& patches = sendSection->getPatches();

     72:         for(typename Section::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
     73:           const Obj<typename Section::section_type>&        section = p_iter->second;
     74:           const typename Section::section_type::chart_type& chart   = section->getChart();

     76:           for(typename Section::section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
     77:             if (sendFiller->hasPoint(*c_iter)) {
     78:               section->updatePoint(*c_iter, sendFiller->restrictPoint(*c_iter));
     79:             }
     80:           }
     81:         }
     82:       }
     83:       template<typename RecvOverlap, typename Sizer, typename Section>
     84:       static void setupReceive(const Obj<RecvOverlap>& recvOverlap, const Obj<Sizer>& recvSizer, const Obj<Section>& recvSection) {
     85:         // Create section
     86:         const Obj<typename RecvOverlap::capSequence> ranks = recvOverlap->cap();

     88:         recvSection->clear();
     89:         for(typename RecvOverlap::capSequence::iterator r_iter = ranks->begin(); r_iter != ranks->end(); ++r_iter) {
     90:           const Obj<typename RecvOverlap::supportSequence>& points  = recvOverlap->support(*r_iter);
     91:           const Obj<typename Section::section_type>&        section = recvSection->getSection(*r_iter);

     93:           // Want to replace this loop with a slice through color
     94:           for(typename RecvOverlap::supportSequence::iterator p_iter = points->begin(); p_iter != points->end(); ++p_iter) {
     95:             const typename dsieve_type::point_type& point = p_iter.color();

     97:             section->setFiberDimension(point, 1);
     98:           }
     99:         }
    100:         recvSection->construct(recvSizer);
    101:         recvSection->allocate();
    102:         recvSection->constructCommunication(Section::RECEIVE);
    103:       }
    104:       template<typename SendOverlap, typename RecvOverlap, typename SizerFiller, typename Filler, typename SendSection, typename RecvSection>
    105:       static void completeSection(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SizerFiller>& sizerFiller, const Filler& filler, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
    106:         typedef typename alloc_type::template rebind<int>::other int_alloc_type;
    107:         typedef typename ALE::Field<SendOverlap, int, ALE::Section<point_type, int, int_alloc_type> > send_sizer_type;
    108:         typedef typename ALE::Field<RecvOverlap, int, ALE::Section<point_type, int, int_alloc_type> > recv_sizer_type;
    109:         typedef typename ALE::Field<SendOverlap, int, ALE::ConstantSection<point_type, int> >         constant_send_sizer;
    110:         typedef typename ALE::Field<RecvOverlap, int, ALE::ConstantSection<point_type, int> >         constant_recv_sizer;
    111:         Obj<send_sizer_type>     sendSizer      = new send_sizer_type(sendSection->comm(), sendSection->debug());
    112:         Obj<recv_sizer_type>     recvSizer      = new recv_sizer_type(recvSection->comm(), sendSizer->getTag(), recvSection->debug());
    113:         Obj<constant_send_sizer> constSendSizer = new constant_send_sizer(sendSection->comm(), sendSection->debug());
    114:         Obj<constant_recv_sizer> constRecvSizer = new constant_recv_sizer(recvSection->comm(), recvSection->debug());
    115:         Obj<topology_type>       sendChart      = completion::createSendTopology(sendOverlap);
    116:         Obj<topology_type>       recvChart      = completion::createRecvTopology(recvOverlap);

    118:         // 1) Create the sizer sections
    119:         constSendSizer->setTopology(sendChart);
    120:         const typename topology_type::sheaf_type& sendRanks = sendChart->getPatches();
    121:         for(typename topology_type::sheaf_type::const_iterator r_iter = sendRanks.begin(); r_iter != sendRanks.end(); ++r_iter) {
    122:           const int rank = r_iter->first;
    123:           const int one  = 1;
    124:           constSendSizer->getSection(rank)->updatePoint(*r_iter->second->base()->begin(), &one);
    125:         }
    126:         constRecvSizer->setTopology(recvChart);
    127:         const typename topology_type::sheaf_type& recvRanks = recvChart->getPatches();
    128:         for(typename topology_type::sheaf_type::const_iterator r_iter = recvRanks.begin(); r_iter != recvRanks.end(); ++r_iter) {
    129:           const int rank = r_iter->first;
    130:           const int one  = 1;
    131:           constRecvSizer->getSection(rank)->updatePoint(*r_iter->second->base()->begin(), &one);
    132:         }
    133:         int_completion::setupSend(sendChart, constSendSizer, sendSizer);
    134:         int_completion::setupReceive(recvOverlap, constRecvSizer, recvSizer);
    135:         // 2) Fill the sizer section and communicate
    136:         int_completion::fillSend(sizerFiller, sendSizer);
    137:         if (sendSizer->debug()) {sendSizer->view("Send Sizer in Completion", MPI_COMM_SELF);}
    138:         sendSizer->startCommunication();
    139:         recvSizer->startCommunication();
    140:         sendSizer->endCommunication();
    141:         recvSizer->endCommunication();
    142:         if (recvSizer->debug()) {recvSizer->view("Receive Sizer in Completion", MPI_COMM_SELF);}
    143:         // No need to update a global section since the receive sizes are all on the interface
    144:         // 3) Create the send and receive sections
    145:         completion::setupSend(sendChart, sendSizer, sendSection);
    146:         completion::setupReceive(recvOverlap, recvSizer, recvSection);
    147:         // 4) Fill up send section and communicate
    148:         completion::fillSend(filler, sendSection);
    149:         if (sendSection->debug()) {sendSection->view("Send Section in Completion", MPI_COMM_SELF);}
    150:         sendSection->startCommunication();
    151:         recvSection->startCommunication();
    152:         sendSection->endCommunication();
    153:         recvSection->endCommunication();
    154:         if (recvSection->debug()) {recvSection->view("Receive Section in Completion", MPI_COMM_SELF);}
    155:       }
    156:     };
    157:   }
    158: }
    159: #endif