Actual source code: SectionCompletion.hh
petsc-3.4.5 2014-06-29
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