source: issm/trunk-jpl/externalpackages/petsc-dev/src/include/petscdmmesh.hh@ 11896

Last change on this file since 11896 was 11896, checked in by habbalf, 13 years ago

petsc-dev : Petsc development code in external packages. Mercurial updates
File size: 82.8 KB
RevLine 
[11896]1#if !defined(__PETSCDMMESH_HH)
2#define __PETSCDMMESH_HH
3
4#include <petscdmmesh.h>
5#include <functional>
6
7using ALE::Obj;
8
9PetscErrorCode DMMeshView_Sieve(const ALE::Obj<PETSC_MESH_TYPE>& mesh, PetscViewer viewer);
10
11#undef __FUNCT__
12#define __FUNCT__ "DMMeshCreateMatrix"
13template<typename Mesh, typename Section>
14PetscErrorCode DMMeshCreateMatrix(const Obj<Mesh>& mesh, const Obj<Section>& section, const MatType mtype, Mat *J, int bs = -1, bool fillMatrix = false)
15{
16 const ALE::Obj<typename Mesh::order_type>& order = mesh->getFactory()->getGlobalOrder(mesh, section->getName(), section);
17 int localSize = order->getLocalSize();
18 int globalSize = order->getGlobalSize();
19 PetscBool isShell, isBlock, isSeqBlock, isMPIBlock, isSymBlock, isSymSeqBlock, isSymMPIBlock, isSymmetric;
20 PetscErrorCode ierr;
21
22 PetscFunctionBegin;
23 ierr = MatCreate(mesh->comm(), J);CHKERRQ(ierr);
24 ierr = MatSetSizes(*J, localSize, localSize, globalSize, globalSize);CHKERRQ(ierr);
25 ierr = MatSetType(*J, mtype);CHKERRQ(ierr);
26 ierr = MatSetFromOptions(*J);CHKERRQ(ierr);
27 ierr = PetscStrcmp(mtype, MATSHELL, &isShell);CHKERRQ(ierr);
28 ierr = PetscStrcmp(mtype, MATBAIJ, &isBlock);CHKERRQ(ierr);
29 ierr = PetscStrcmp(mtype, MATSEQBAIJ, &isSeqBlock);CHKERRQ(ierr);
30 ierr = PetscStrcmp(mtype, MATMPIBAIJ, &isMPIBlock);CHKERRQ(ierr);
31 ierr = PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);CHKERRQ(ierr);
32 ierr = PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);CHKERRQ(ierr);
33 ierr = PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);CHKERRQ(ierr);
34 // Check for symmetric storage
35 isSymmetric = (PetscBool) (isSymBlock || isSymSeqBlock || isSymMPIBlock);
36 if (!isShell) {
37 PetscInt *dnz, *onz, bsLocal;
38
39 if (bs < 0) {
40 if (isBlock || isSeqBlock || isMPIBlock || isSymBlock || isSymSeqBlock || isSymMPIBlock) {
41 const typename Section::chart_type& chart = section->getChart();
42
43 for(typename Section::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
44 if (section->getFiberDimension(*c_iter)) {
45 bs = section->getFiberDimension(*c_iter);
46 break;
47 }
48 }
49 } else {
50 bs = 1;
51 }
52 // Must have same blocksize on all procs (some might have no points)
53 bsLocal = bs;
54 ierr = MPI_Allreduce(&bsLocal, &bs, 1, MPIU_INT, MPI_MAX, mesh->comm());CHKERRQ(ierr);
55 }
56 ierr = PetscMalloc2(localSize/bs, PetscInt, &dnz, localSize/bs, PetscInt, &onz);CHKERRQ(ierr);
57#ifdef USE_NEW_OVERLAP
58 ierr = preallocateOperatorNewOverlap(mesh, bs, section->getAtlas(), order, dnz, onz, isSymmetric, *J, fillMatrix);CHKERRQ(ierr);
59#else
60 ierr = preallocateOperatorNew(mesh, bs, section->getAtlas(), order, dnz, onz, isSymmetric, *J, fillMatrix);CHKERRQ(ierr);
61#endif
62 ierr = PetscFree2(dnz, onz);CHKERRQ(ierr);
63 if (isSymmetric) {
64 ierr = MatSetOption(*J, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);CHKERRQ(ierr);
65 }
66 }
67 PetscFunctionReturn(0);
68}
69
70#undef __FUNCT__
71#define __FUNCT__ "DMMeshCreateGlobalScatter"
72template<typename Mesh, typename Section>
73PetscErrorCode DMMeshCreateGlobalScatter(const ALE::Obj<Mesh>& m, const ALE::Obj<Section>& s, VecScatter *scatter)
74{
75 const ALE::Obj<typename Mesh::order_type>& globalOrder = m->getFactory()->getGlobalOrder(m, s->getName(), s);
76 PetscErrorCode ierr;
77
78 PetscFunctionBegin;
79 ierr = DMMeshCreateGlobalScatter(m, s, globalOrder, false, scatter);CHKERRQ(ierr);
80 PetscFunctionReturn(0);
81}
82
83#undef __FUNCT__
84#define __FUNCT__ "DMMeshCreateGlobalScatter"
85template<typename Mesh, typename Section>
86PetscErrorCode DMMeshCreateGlobalScatter(const ALE::Obj<Mesh>& m, const ALE::Obj<Section>& s, const ALE::Obj<typename Mesh::label_type>& label, VecScatter *scatter)
87{
88 const ALE::Obj<typename Mesh::order_type>& globalOrder = m->getFactory()->getGlobalOrder(m, s->getName(), s, -1, label);
89 PetscErrorCode ierr;
90
91 PetscFunctionBegin;
92 ierr = DMMeshCreateGlobalScatter(m, s, globalOrder, false, scatter);CHKERRQ(ierr);
93 PetscFunctionReturn(0);
94}
95
96#undef __FUNCT__
97#define __FUNCT__ "DMMeshCreateGlobalScatter"
98template<typename Mesh, typename Section>
99PetscErrorCode DMMeshCreateGlobalScatter(const ALE::Obj<Mesh>& m, const std::string& name, const typename Section::chart_type& points, const ALE::Obj<Section>& s, VecScatter *scatter)
100{
101 const ALE::Obj<typename Mesh::order_type>& globalOrder = m->getFactory()->getGlobalOrder(m, name, points, s);
102 PetscErrorCode ierr;
103
104 PetscFunctionBegin;
105 ierr = DMMeshCreateGlobalScatter(m, s, globalOrder, false, scatter);CHKERRQ(ierr);
106 PetscFunctionReturn(0);
107}
108
109#undef __FUNCT__
110#define __FUNCT__ "DMMeshCreateGlobalScatter"
111template<typename Mesh, typename Section>
112PetscErrorCode DMMeshCreateGlobalScatter(const ALE::Obj<Mesh>& m, const ALE::Obj<Section>& s, const ALE::Obj<typename Mesh::order_type>& globalOrder, bool includeConstraints, VecScatter *scatter)
113{
114 typedef typename Mesh::real_section_type::index_type index_type;
115 PetscErrorCode ierr;
116
117 PetscFunctionBegin;
118 ierr = PetscLogEventBegin(DMMesh_GetGlobalScatter,0,0,0,0);CHKERRQ(ierr);
119 const typename Mesh::order_type::chart_type& chart = globalOrder->getChart();
120 int *localIndices, *globalIndices;
121 int localSize = globalOrder->getLocalSize();
122 int overlapSize = -1;
123 int localIndx = 0, globalIndx = 0;
124 Vec globalVec, localVec;
125 IS localIS, globalIS;
126
127 ierr = VecCreate(m->comm(), &globalVec);CHKERRQ(ierr);
128 ierr = VecSetSizes(globalVec, localSize, PETSC_DETERMINE);CHKERRQ(ierr);
129 ierr = VecSetFromOptions(globalVec);CHKERRQ(ierr);
130
131 if (includeConstraints) {
132 overlapSize = s->sizeWithBC();
133 ierr = PetscMalloc(overlapSize*sizeof(int), &localIndices);CHKERRQ(ierr);
134 ierr = PetscMalloc(overlapSize*sizeof(int), &globalIndices);CHKERRQ(ierr);
135 } else {
136 overlapSize = s->size();
137 ierr = PetscMalloc(overlapSize*sizeof(int), &localIndices);CHKERRQ(ierr);
138 ierr = PetscMalloc(overlapSize*sizeof(int), &globalIndices);CHKERRQ(ierr);
139 } // if/else
140
141 // Loop over all local points
142 for(typename Mesh::order_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
143 // Map local indices to global indices
144 if (includeConstraints) {
145 s->getIndicesRaw(*p_iter, localIndices, &localIndx, 0);
146 s->getIndicesRaw(*p_iter, globalOrder, globalIndices, &globalIndx, 0);
147 } else {
148 s->getIndices(*p_iter, localIndices, &localIndx, 0, true, true);
149 s->getIndices(*p_iter, globalOrder, globalIndices, &globalIndx, 0, true, false);
150 }
151 //numConstraints += s->getConstraintDimension(*p_iter);
152 }
153 // Local arrays also have constraints, which are not mapped
154 if (localIndx > overlapSize) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Invalid number of local indices %d, should not be greater than %d", localIndx, overlapSize);
155 if (globalIndx > overlapSize) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Invalid number of global indices %d, should not be greater than %d", globalIndx, overlapSize);
156 if (globalIndx != localIndx) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Mismatched number of global indices %d, and local indices %d", globalIndx, localIndx);
157 if (m->debug()) {
158 globalOrder->view("Global Order");
159 for(int i = 0; i < globalIndx; ++i) {
160 printf("[%d] localIndex[%d]: %d globalIndex[%d]: %d\n", m->commRank(), i, localIndices[i], i, globalIndices[i]);
161 }
162 }
163 ierr = ISCreateGeneral(PETSC_COMM_SELF, localIndx, localIndices,PETSC_OWN_POINTER, &localIS);CHKERRQ(ierr);
164 ierr = ISCreateGeneral(PETSC_COMM_SELF, globalIndx, globalIndices,PETSC_OWN_POINTER, &globalIS);CHKERRQ(ierr);
165 // Can remove this when I test it with NULL
166#ifdef PETSC_USE_COMPLEX
167 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF, s->getStorageSize(), PETSC_NULL, &localVec);CHKERRQ(ierr);
168#else
169 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF, s->getStorageSize(), s->restrictSpace(), &localVec);CHKERRQ(ierr);
170#endif
171 ierr = VecScatterCreate(localVec, localIS, globalVec, globalIS, scatter);CHKERRQ(ierr);
172 ierr = ISDestroy(&globalIS);CHKERRQ(ierr);
173 ierr = ISDestroy(&localIS);CHKERRQ(ierr);
174 ierr = VecDestroy(&localVec);CHKERRQ(ierr);
175 ierr = VecDestroy(&globalVec);CHKERRQ(ierr);
176 ierr = PetscLogEventEnd(DMMesh_GetGlobalScatter,0,0,0,0);CHKERRQ(ierr);
177 PetscFunctionReturn(0);
178}
179
180template<typename Mesh, typename Section>
181void createOperator(const ALE::Obj<Mesh>& mesh, const ALE::Obj<Section>& s, const ALE::Obj<Mesh>& op) {
182 typedef ALE::SieveAlg<Mesh> sieve_alg_type;
183 typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
184 const typename Section::chart_type& chart = s->getChart();
185
186 // Create local operator
187 // We do not decorate arrows yet
188 for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
189 const Obj<typename sieve_alg_type::supportArray>& star = sieve_alg_type::star(mesh, *p_iter);
190
191 for(typename sieve_alg_type::supportArray::const_iterator s_iter = star->begin(); s_iter != star->end(); ++s_iter) {
192 const Obj<typename sieve_alg_type::coneArray>& closure = sieve_alg_type::closure(mesh, *s_iter);
193
194 for(typename sieve_alg_type::coneArray::const_iterator c_iter = closure->begin(); c_iter != closure->end(); ++c_iter) {
195 op->getSieve()->addCone(*c_iter, *p_iter);
196 }
197 }
198 }
199 op->view("Local operator");
200 // Construct overlap
201 Obj<FlexMesh::send_overlap_type> sendOverlap = mesh->getSendOverlap();
202 Obj<FlexMesh::recv_overlap_type> recvOverlap = mesh->getRecvOverlap();
203 FlexMesh::renumbering_type& renumbering = mesh->getRenumbering();
204
205 sendOverlap->view("Mesh send overlap");
206 recvOverlap->view("Mesh recv overlap");
207 // Complete operator
208 typedef ALE::DistributionNew<FlexMesh>::cones_type ConeOverlap;
209
210 ALE::Obj<ConeOverlap> overlapCones = ALE::DistributionNew<FlexMesh>::completeCones(op->getSieve(), op->getSieve(), renumbering, sendOverlap, recvOverlap);
211 op->view("Completed operator");
212 // Update renumbering and overlap
213 overlapCones->view("Overlap cones");
214 Obj<FlexMesh::send_overlap_type> opSendOverlap = op->getSendOverlap();
215 Obj<FlexMesh::recv_overlap_type> opRecvOverlap = op->getRecvOverlap();
216 FlexMesh::renumbering_type& opRenumbering = op->getRenumbering();
217 const typename ConeOverlap::chart_type& overlapChart = overlapCones->getChart();
218 int p = renumbering.size();
219
220 opRenumbering = renumbering;
221 for(typename ConeOverlap::chart_type::const_iterator p_iter = overlapChart.begin(); p_iter != overlapChart.end(); ++p_iter) {
222 if (opRenumbering.find(p_iter->second) == opRenumbering.end()) {
223 opRenumbering[p_iter->second] = p++;
224 }
225 }
226 ALE::SetFromMap<FlexMesh::renumbering_type> opGlobalPoints(opRenumbering);
227
228 ALE::OverlapBuilder<>::constructOverlap(opGlobalPoints, opRenumbering, opSendOverlap, opRecvOverlap);
229 sendOverlap->view("Operator send overlap");
230 recvOverlap->view("Operator recv overlap");
231 // Create global order
232 Obj<FlexMesh::order_type> globalOrder = new FlexMesh::order_type(op->comm(), op->debug());
233
234 op->getFactory()->constructLocalOrder(globalOrder, opSendOverlap, opGlobalPoints, s);
235 op->getFactory()->calculateOffsets(globalOrder);
236 op->getFactory()->updateOrder(globalOrder, opGlobalPoints);
237 op->getFactory()->completeOrder(globalOrder, opSendOverlap, opRecvOverlap);
238 globalOrder->view("Operator global order");
239 // Create dnz/onz or CSR
240};
241
242template<typename Atlas>
243class AdjVisitor {
244public:
245 typedef typename ALE::Mesh<PetscInt,PetscScalar>::point_type point_type;
246protected:
247 Atlas& atlas;
248 ALE::Mesh<PetscInt,PetscScalar>::sieve_type& adjGraph;
249 point_type p;
250public:
251 AdjVisitor(Atlas& atlas, ALE::Mesh<PetscInt,PetscScalar>::sieve_type& adjGraph) : atlas(atlas), adjGraph(adjGraph) {};
252 void visitPoint(const point_type& point) {
253 if (atlas.restrictPoint(point)[0].prefix) {
254 adjGraph.addCone(point, p);
255 }
256 };
257 template<typename Arrow>
258 void visitArrow(const Arrow&) {};
259public:
260 void setRoot(const point_type& point) {this->p = point;};
261};
262
263#undef __FUNCT__
264#define __FUNCT__ "createLocalAdjacencyGraph"
265template<typename Mesh, typename Atlas>
266PetscErrorCode createLocalAdjacencyGraph(const ALE::Obj<Mesh>& mesh, const ALE::Obj<Atlas>& atlas, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::sieve_type>& adjGraph)
267{
268 typedef typename ALE::ISieveVisitor::TransitiveClosureVisitor<typename Mesh::sieve_type, AdjVisitor<Atlas> > ClosureVisitor;
269 typedef typename ALE::ISieveVisitor::ConeVisitor<typename Mesh::sieve_type, ClosureVisitor> ConeVisitor;
270 typedef typename ALE::ISieveVisitor::TransitiveClosureVisitor<typename Mesh::sieve_type, ConeVisitor> StarVisitor;
271 AdjVisitor<Atlas> adjV(*atlas, *adjGraph);
272 ClosureVisitor closureV(*mesh->getSieve(), adjV);
273 ConeVisitor coneV(*mesh->getSieve(), closureV);
274 StarVisitor starV(*mesh->getSieve(), coneV);
275 /* In general, we need to get FIAT info that attaches dual basis vectors to sieve points */
276 const typename Atlas::chart_type& chart = atlas->getChart();
277
278 PetscFunctionBegin;
279 starV.setIsCone(false);
280 for(typename Atlas::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
281 adjV.setRoot(*c_iter);
282 mesh->getSieve()->support(*c_iter, starV);
283 closureV.clear();
284 starV.clear();
285 }
286 PetscFunctionReturn(0);
287}
288
289#undef __FUNCT__
290#define __FUNCT__ "createLocalAdjacencyGraphI"
291template<typename Mesh, typename Atlas>
292PetscErrorCode createLocalAdjacencyGraphI(const ALE::Obj<Mesh>& mesh, const ALE::Obj<Atlas>& atlas, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::sieve_type>& adjGraph)
293{
294 typedef typename ALE::ISieveVisitor::TransitiveClosureVisitor<typename Mesh::sieve_type, AdjVisitor<Atlas> > ClosureVisitor;
295 typedef typename ALE::ISieveVisitor::ConeVisitor<typename Mesh::sieve_type, ClosureVisitor> ConeVisitor;
296 typedef typename ALE::ISieveVisitor::TransitiveClosureVisitor<typename Mesh::sieve_type, ConeVisitor> StarVisitor;
297 AdjVisitor<Atlas> adjV(*atlas, *adjGraph);
298 ClosureVisitor closureV(*mesh->getSieve(), adjV);
299 ConeVisitor coneV(*mesh->getSieve(), closureV);
300 StarVisitor starV(*mesh->getSieve(), coneV);
301 /* In general, we need to get FIAT info that attaches dual basis vectors to sieve points */
302 const typename Atlas::chart_type& chart = atlas->getChart();
303
304 PetscFunctionBegin;
305 // Changes for ISieve
306 // 1) Add AdjSizeVisitor to set cone sizes
307 // 2) Add new symmetrizeSizes() to ISieve to calculate support sizes
308 // 3) Allocate adjGraph
309 // 4) Change AdjVisitor to stack up cone rather than calling addPoint()
310 // 5) Get points and call setCone() each time
311 starV.setIsCone(false);
312 for(typename Atlas::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
313 adjV.setRoot(*c_iter);
314 mesh->getSieve()->support(*c_iter, starV);
315 closureV.clear();
316 starV.clear();
317 }
318 PetscFunctionReturn(0);
319}
320
321#undef __FUNCT__
322#define __FUNCT__ "createLocalAdjacencyGraph"
323template<typename Atlas>
324PetscErrorCode createLocalAdjacencyGraph(const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar> >& mesh, const ALE::Obj<Atlas>& atlas, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::sieve_type>& adjGraph)
325{
326 typedef ALE::SieveAlg<ALE::Mesh<PetscInt,PetscScalar> > sieve_alg_type;
327 /* In general, we need to get FIAT info that attaches dual basis vectors to sieve points */
328 const typename Atlas::chart_type& chart = atlas->getChart();
329
330 PetscFunctionBegin;
331 for(typename Atlas::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
332 const Obj<typename sieve_alg_type::supportArray>& star = sieve_alg_type::star(mesh, *c_iter);
333
334 for(typename sieve_alg_type::supportArray::const_iterator s_iter = star->begin(); s_iter != star->end(); ++s_iter) {
335 const Obj<typename sieve_alg_type::coneArray>& closure = sieve_alg_type::closure(mesh, *s_iter);
336
337 for(typename sieve_alg_type::coneArray::const_iterator cl_iter = closure->begin(); cl_iter != closure->end(); ++cl_iter) {
338 adjGraph->addCone(*cl_iter, *c_iter);
339 }
340 }
341 }
342 PetscFunctionReturn(0);
343}
344
345template<typename Arrow>
346struct SelectSource : public std::unary_function<Arrow, typename Arrow::source_type>
347{
348 typename Arrow::source_type& operator()(Arrow& x) const {return x.source;}
349 const typename Arrow::source_type& operator()(const Arrow& x) const {return x.source;}
350};
351
352template<class Pair>
353struct Select1st : public std::unary_function<Pair, typename Pair::first_type>
354{
355 typename Pair::first_type& operator()(Pair& x) const {return x.first;}
356 const typename Pair::first_type& operator()(const Pair& x) const {return x.first;}
357};
358
359template<typename Mesh, typename Order>
360PetscErrorCode globalizeLocalAdjacencyGraph(const ALE::Obj<Mesh>& mesh, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::sieve_type>& adjGraph, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::send_overlap_type>& sendOverlap, const ALE::Obj<Order>& globalOrder)
361{
362 typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
363 const int debug = mesh->debug();
364
365 PetscFunctionBegin;
366 ALE::PointFactory<typename Mesh::point_type>& pointFactory = ALE::PointFactory<FlexMesh::point_type>::singleton(mesh->comm(), mesh->getSieve()->getChart().max(), mesh->debug());
367 // Check for points not in sendOverlap
368 std::set<typename Mesh::point_type> interiorPoints;
369 std::set<typename Mesh::point_type> overlapSources;
370 std::set<typename Mesh::sieve_type::arrow_type> overlapArrows;
371 const Obj<FlexMesh::sieve_type::traits::capSequence>& columns = adjGraph->cap();
372
373 for(FlexMesh::sieve_type::traits::capSequence::iterator p_iter = columns->begin(); p_iter != columns->end(); ++p_iter) {
374 // This optimization does not work because empty points are sent anyway
375 //if (!sendOverlap->support(*p_iter)->size() && globalOrder->restrictPoint(*p_iter)[0].index) {
376 if (!sendOverlap->support(*p_iter)->size()) {
377 interiorPoints.insert(*p_iter);
378 } else {
379 // If a point p is in the overlap for process i and an adjacent point q is in the overlap for process j then:
380 // Replace (q, p) with (q', p)
381 // Notice I can reverse the arrow because the graph is initially symmetric
382 if (debug) {std::cout << "["<<globalOrder->commRank()<<"] Checking overlap point " << *p_iter << " for overlap neighbors" << std::endl;}
383 const Obj<typename FlexMesh::sieve_type::supportSequence>& neighbors = adjGraph->support(*p_iter);
384 const typename FlexMesh::sieve_type::supportSequence::iterator nEnd = neighbors->end();
385
386 for(typename FlexMesh::sieve_type::supportSequence::iterator n_iter = neighbors->begin(); n_iter != nEnd; ++n_iter) {
387 if (sendOverlap->support(*n_iter)->size()) {
388 if (debug) {std::cout << "["<<globalOrder->commRank()<<"] Found overlap neighbor " << *n_iter << std::endl;}
389 const Obj<typename FlexMesh::send_overlap_type::supportSequence>& ranks = sendOverlap->support(*p_iter);
390 const typename FlexMesh::send_overlap_type::supportSequence::iterator rEnd = ranks->end();
391 bool equal = true;
392
393 for(typename FlexMesh::send_overlap_type::supportSequence::iterator r_iter = ranks->begin(); r_iter != rEnd; ++r_iter) {
394 const Obj<typename FlexMesh::send_overlap_type::supportSequence>& nRanks = sendOverlap->support(*n_iter);
395 const typename FlexMesh::send_overlap_type::supportSequence::iterator nrEnd = nRanks->end();
396 bool found = false;
397
398 if (debug) {std::cout << "["<<globalOrder->commRank()<<"] Checking overlap rank " << *r_iter << std::endl;}
399 for(typename FlexMesh::send_overlap_type::supportSequence::iterator nr_iter = nRanks->begin(); nr_iter != nrEnd; ++nr_iter) {
400 if (debug) {std::cout << "["<<globalOrder->commRank()<<"] Checking neighbor overlap rank " << *nr_iter << std::endl;}
401 if (*nr_iter == *r_iter) {
402 found = true;
403 break;
404 }
405 }
406 if (!found) {
407 equal = false;
408 break;
409 }
410 }
411 if (!equal) {
412 if (debug) {std::cout << "["<<globalOrder->commRank()<<"] Unequal rank sets, replacing arrow " << *n_iter <<" --> "<<*p_iter << std::endl;}
413 overlapArrows.insert(typename Mesh::sieve_type::arrow_type(*n_iter, *p_iter));
414 } else {
415 if (debug) {std::cout << "["<<globalOrder->commRank()<<"] Equal rank sets, doing nothing for arrow " << *n_iter <<" --> "<<*p_iter << std::endl;}
416 }
417 }
418 }
419 }
420 }
421 // Renumber those points
422 pointFactory.clear();
423 pointFactory.setMax(mesh->getSieve()->getChart().max());
424 pointFactory.renumberPoints(interiorPoints.begin(), interiorPoints.end());
425 //pointFactory.renumberPoints(overlapArrows.begin(), overlapArrows.end(), SelectSource<typename Mesh::sieve_type::arrow_type>());
426 // They should use a key extractor: overlapSources.insert(overlapArrows.begin(), overlapArrows.end(), SelectSource<typename Mesh::sieve_type::arrow_type>());
427 for(typename std::set<typename Mesh::sieve_type::arrow_type>::const_iterator a_iter = overlapArrows.begin(); a_iter != overlapArrows.end(); ++a_iter) {
428 overlapSources.insert(a_iter->source);
429 }
430 pointFactory.renumberPoints(overlapSources.begin(), overlapSources.end());
431 typename ALE::PointFactory<typename Mesh::point_type>::renumbering_type& renumbering = pointFactory.getRenumbering();
432 typename ALE::PointFactory<typename Mesh::point_type>::renumbering_type& invRenumbering = pointFactory.getInvRenumbering();
433 // Replace points in local adjacency graph
434 const Obj<FlexMesh::sieve_type::traits::baseSequence>& base = adjGraph->base();
435 ALE::Obj<std::vector<typename Mesh::point_type> > newCone = new std::vector<typename Mesh::point_type>();
436
437 for(FlexMesh::sieve_type::traits::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
438 const ALE::Obj<FlexMesh::sieve_type::coneSequence>& cone = adjGraph->cone(*b_iter);
439 const FlexMesh::sieve_type::coneSequence::iterator cEnd = cone->end();
440 bool replace = false;
441
442 for(FlexMesh::sieve_type::coneSequence::iterator c_iter = cone->begin(); c_iter != cEnd; ++c_iter) {
443 if (interiorPoints.find(*c_iter) != interiorPoints.end()) {
444 newCone->push_back(invRenumbering[*c_iter]);
445 replace = true;
446 } else {
447 newCone->push_back(*c_iter);
448 }
449 }
450 if (interiorPoints.find(*b_iter) != interiorPoints.end()) {
451 adjGraph->clearCone(*b_iter);
452 adjGraph->setCone(newCone, invRenumbering[*b_iter]);
453 if (debug) {std::cout << "["<<globalOrder->commRank()<<"]: Replacing cone for " << *b_iter << "("<<invRenumbering[*b_iter]<<") with" << std::endl;}
454 } else if (replace) {
455 adjGraph->clearCone(*b_iter);
456 adjGraph->setCone(newCone, *b_iter);
457 if (debug) {std::cout << "["<<globalOrder->commRank()<<"]: Replacing cone for " << *b_iter << " with" << std::endl;}
458 }
459 if (debug && ((interiorPoints.find(*b_iter) != interiorPoints.end()) || replace)) {
460 for(typename std::vector<typename Mesh::point_type>::const_iterator c_iter = newCone->begin(); c_iter != newCone->end(); ++c_iter) {
461 std::cout << "["<<globalOrder->commRank()<<"]: point " << *c_iter << std::endl;
462 }
463 }
464 newCone->clear();
465 }
466 // Replace arrows
467 for(typename std::set<typename Mesh::sieve_type::arrow_type>::const_iterator a_iter = overlapArrows.begin(); a_iter != overlapArrows.end(); ++a_iter) {
468 if (debug) {std::cout << "["<<globalOrder->commRank()<<"]: Replacing " << a_iter->source<<" --> "<<a_iter->target << " with " << invRenumbering[a_iter->source]<<" --> "<<a_iter->target << std::endl;}
469 adjGraph->removeArrow(a_iter->source, a_iter->target);
470 adjGraph->addArrow(invRenumbering[a_iter->source], a_iter->target);
471 }
472 // Add new points into ordering
473#if 1
474 for(typename ALE::PointFactory<typename Mesh::point_type>::renumbering_type::const_iterator p_iter = renumbering.begin(); p_iter != renumbering.end(); ++p_iter) {
475 if (debug) {std::cout << "["<<globalOrder->commRank()<<"]: Updating " << p_iter->first << " to " << globalOrder->restrictPoint(p_iter->second)[0] << std::endl;}
476 globalOrder->addPoint(p_iter->first);
477 globalOrder->updatePoint(p_iter->first, globalOrder->restrictPoint(p_iter->second));
478 }
479#else
480 globalOrder->reallocatePoint(renumbering.begin(), renumbering.end(), Select1st<typename ALE::PointFactory<typename Mesh::point_type>::renumbering_type::const_iterator::value_type>());
481 for(typename ALE::PointFactory<typename Mesh::point_type>::renumbering_type::const_iterator p_iter = renumbering.begin(); p_iter != renumbering.end(); ++p_iter) {
482 if (debug) {std::cout << "["<<globalOrder->commRank()<<"]: Updating " << p_iter->first << " to " << globalOrder->restrictPoint(p_iter->second)[0] << std::endl;}
483 globalOrder->updatePoint(p_iter->first, globalOrder->restrictPoint(p_iter->second));
484 }
485#endif
486 PetscFunctionReturn(0);
487}
488
489template<typename Order>
490PetscErrorCode globalizeLocalAdjacencyGraph(const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar> >& mesh, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::sieve_type>& adjGraph, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::send_overlap_type>& sendOverlap, const ALE::Obj<Order>& globalOrder)
491{
492 PetscFunctionBegin;
493 PetscFunctionReturn(0);
494}
495
496template<typename Mesh, typename Order>
497PetscErrorCode localizeLocalAdjacencyGraph(const ALE::Obj<Mesh>& mesh, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::sieve_type>& adjGraph, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::send_overlap_type>& sendOverlap, const ALE::Obj<Order>& globalOrder)
498{
499 PetscFunctionBegin;
500 ALE::PointFactory<typename Mesh::point_type>& pointFactory = ALE::PointFactory<ALE::Mesh<PetscInt,PetscScalar>::point_type>::singleton(mesh->comm(), mesh->getSieve()->getChart().max(), mesh->debug());
501 typename ALE::PointFactory<typename Mesh::point_type>::renumbering_type& renumbering = pointFactory.getRenumbering();
502 // Replace points in local adjacency graph
503 const Obj<ALE::Mesh<PetscInt,PetscScalar>::sieve_type::traits::baseSequence>& base = adjGraph->base();
504
505 for(ALE::Mesh<PetscInt,PetscScalar>::sieve_type::traits::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
506 // Replace globalized points
507 if (renumbering.find(*b_iter) != renumbering.end()) {
508 adjGraph->clearCone(renumbering[*b_iter]);
509 adjGraph->setCone(adjGraph->cone(*b_iter), renumbering[*b_iter]);
510 adjGraph->clearCone(*b_iter);
511 }
512 }
513 PetscFunctionReturn(0);
514}
515
516template<typename Order>
517PetscErrorCode localizeLocalAdjacencyGraph(const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar> >& mesh, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::sieve_type>& adjGraph, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::send_overlap_type>& sendOverlap, const ALE::Obj<Order>& globalOrder)
518{
519 PetscFunctionBegin;
520 PetscFunctionReturn(0);
521}
522
523template<typename Mesh, typename Order>
524PetscErrorCode renumberLocalAdjacencyGraph(const ALE::Obj<Mesh>& mesh, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::sieve_type>& adjGraph, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::send_overlap_type>& sendOverlap, const ALE::Obj<Order>& globalOrder)
525{
526 typedef typename ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
527 ALE::Obj<std::set<typename Mesh::point_type> > newCone = new std::set<typename Mesh::point_type>();
528 const int debug = mesh->debug();
529
530 PetscFunctionBegin;
531 ALE::PointFactory<typename Mesh::point_type>& pointFactory = ALE::PointFactory<typename FlexMesh::point_type>::singleton(mesh->comm(), mesh->getSieve()->getChart().max(), debug);
532 pointFactory.constructRemoteRenumbering();
533 typename ALE::PointFactory<typename Mesh::point_type>::renumbering_type& renumbering = pointFactory.getRenumbering();
534 typename ALE::PointFactory<typename Mesh::point_type>::remote_renumbering_type& remoteRenumbering = pointFactory.getRemoteRenumbering();
535 // Replace points in local adjacency graph
536 const Obj<typename FlexMesh::sieve_type::traits::baseSequence>& base = adjGraph->base();
537
538 for(FlexMesh::sieve_type::traits::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
539 // Loop over cone checking for remote points that shadow local points
540 const Obj<FlexMesh::sieve_type::traits::coneSequence>& cone = adjGraph->cone(*b_iter);
541 const FlexMesh::sieve_type::traits::coneSequence::iterator cEnd = cone->end();
542 bool replace = false;
543
544 if (debug) {std::cout <<"["<<adjGraph->commRank()<<"]: Checking base point " << *b_iter << std::endl;}
545 for(FlexMesh::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cEnd; ++c_iter) {
546 bool useOldPoint = true;
547
548 if (debug) {std::cout <<"["<<adjGraph->commRank()<<"]: cone point " << *c_iter;}
549 if (!globalOrder->isLocal(*c_iter)) {
550 if (debug) {std::cout << " is nonlocal" << std::endl;}
551 for(int p = 0; p < adjGraph->commSize(); ++p) {
552 if (p == adjGraph->commRank()) continue;
553 if (remoteRenumbering[p].find(*c_iter) != remoteRenumbering[p].end()) {
554 if (debug) {std::cout <<"["<<adjGraph->commRank()<<"]: found " << *c_iter << " on process " << p << " as point " << remoteRenumbering[p][*c_iter];}
555 const Obj<FlexMesh::send_overlap_type::traits::coneSequence>& localPoint = sendOverlap->cone(p, remoteRenumbering[p][*c_iter]);
556
557 if (localPoint->size()) {
558 if (debug) {std::cout << " with local match " << *localPoint->begin() << std::endl;}
559 newCone->insert(*localPoint->begin());
560 replace = true;
561 useOldPoint = false;
562 break;
563 } else {
564 if (debug) {std::cout << " but not in sendOverlap" << std::endl;}
565 }
566 }
567 }
568 } else {
569 if (debug) {std::cout << " is local" << std::endl;}
570 if (renumbering.find(*c_iter) != renumbering.end()) {
571 if (debug) {std::cout <<"["<<adjGraph->commRank()<<"]: found " << *c_iter << " locally as point " << renumbering[*c_iter];}
572 newCone->insert(renumbering[*c_iter]);
573 replace = true;
574 useOldPoint = false;
575 }
576 }
577 if (useOldPoint) {
578 newCone->insert(*c_iter);
579 }
580 }
581 if (replace) {
582 if (debug > 1) {
583 std::cout <<"["<<adjGraph->commRank()<<"]: Replacing cone for " << *b_iter << " due to shadowed points from" << std::endl;
584 const Obj<FlexMesh::sieve_type::traits::coneSequence>& cone = adjGraph->cone(*b_iter);
585 const FlexMesh::sieve_type::traits::coneSequence::iterator cEnd = cone->end();
586 for(typename FlexMesh::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cEnd; ++c_iter) {
587 std::cout <<"["<<adjGraph->commRank()<<"]: point " << *c_iter << std::endl;
588 }
589 std::cout <<"["<<adjGraph->commRank()<<"]: to" << std::endl;
590 for(typename std::set<typename Mesh::point_type>::const_iterator c_iter = newCone->begin(); c_iter != newCone->end(); ++c_iter) {
591 std::cout <<"["<<adjGraph->commRank()<<"]: point " << *c_iter << std::endl;
592 }
593 }
594 adjGraph->setCone(newCone, *b_iter);
595 newCone->clear();
596 }
597 }
598 PetscFunctionReturn(0);
599}
600
601template<typename Order>
602PetscErrorCode renumberLocalAdjacencyGraph(const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar> >& mesh, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::sieve_type>& adjGraph, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::send_overlap_type>& sendOverlap, const ALE::Obj<Order>& globalOrder)
603{
604 PetscFunctionBegin;
605 PetscFunctionReturn(0);
606}
607
608#undef __FUNCT__
609#define __FUNCT__ "preallocateOperator"
610/* Problem:
611 We want to allocate an operator. This means we want to know the ordering of all unknowns in the sparsity pattern.
612 The preexisting overlap will not contain information about all unknowns (columns) in the operator.
613
614 Solution:
615 Construct the local sparsity pattern, using globally consistent names for new interior points. Cone complete this
616 pattern, which gives an augmented overlap structure. Insert offsets for the new, global point ids in the existing
617 order, and then complete the order. This argues for a recreation of the order, rather than use of an existing
618 order.
619
620 Problem:
621 Figure out sparsity pattern of the operator, when we have already locally numbered all points. The overlap can
622 establish common names for shared points, but not for interior points.
623
624 Solution:
625 Create a shared resource that bestows globally consistent names.
626*/
627template<typename Mesh, typename Atlas>
628PetscErrorCode preallocateOperator(const ALE::Obj<Mesh>& mesh, const int bs, const ALE::Obj<Atlas>& atlas, const ALE::Obj<typename Mesh::order_type>& globalOrder, PetscInt dnz[], PetscInt onz[], Mat A)
629{
630 MPI_Comm comm = mesh->comm();
631 const int rank = mesh->commRank();
632 const int debug = mesh->debug();
633 const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar> > adjBundle = new ALE::Mesh<PetscInt,PetscScalar>(comm, debug);
634 const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::sieve_type> adjGraph = new ALE::Mesh<PetscInt,PetscScalar>::sieve_type(comm, debug);
635 PetscInt numLocalRows, firstRow;
636 PetscErrorCode ierr;
637
638 PetscFunctionBegin;
639 adjBundle->setSieve(adjGraph);
640 numLocalRows = globalOrder->getLocalSize();
641 firstRow = globalOrder->getGlobalOffsets()[rank];
642 ierr = createLocalAdjacencyGraph(mesh, atlas, adjGraph);CHKERRQ(ierr);
643 /* Distribute adjacency graph */
644 adjBundle->constructOverlap();
645 typedef typename Mesh::sieve_type::point_type point_type;
646 typedef typename Mesh::send_overlap_type send_overlap_type;
647 typedef typename Mesh::recv_overlap_type recv_overlap_type;
648 typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, point_type> > send_section_type;
649 typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, point_type> > recv_section_type;
650 const Obj<send_overlap_type>& vertexSendOverlap = mesh->getSendOverlap();
651 const Obj<recv_overlap_type>& vertexRecvOverlap = mesh->getRecvOverlap();
652 const Obj<send_overlap_type> nbrSendOverlap = new send_overlap_type(comm, debug);
653 const Obj<recv_overlap_type> nbrRecvOverlap = new recv_overlap_type(comm, debug);
654 const Obj<send_section_type> sendSection = new send_section_type(comm, debug);
655 const Obj<recv_section_type> recvSection = new recv_section_type(comm, sendSection->getTag(), debug);
656
657 if (mesh->commSize() > 1) {
658 if (debug > 1) adjGraph->view("Original adjacency graph");
659 ierr = globalizeLocalAdjacencyGraph(mesh, adjGraph, vertexSendOverlap, globalOrder);
660 if (debug > 1) adjGraph->view("Globalized adjacency graph");
661 ALE::Distribution<ALE::Mesh<PetscInt,PetscScalar> >::coneCompletion(vertexSendOverlap, vertexRecvOverlap, adjBundle, sendSection, recvSection);
662 if (debug > 1) adjGraph->view("Completed adjacency graph");
663 ierr = localizeLocalAdjacencyGraph(mesh, adjGraph, vertexSendOverlap, globalOrder);
664 if (debug > 1) adjGraph->view("Localized adjacency graph");
665 /* Distribute indices for new points */
666 ALE::Distribution<ALE::Mesh<PetscInt,PetscScalar> >::updateOverlap(vertexSendOverlap, vertexRecvOverlap, sendSection, recvSection, nbrSendOverlap, nbrRecvOverlap);
667 mesh->getFactory()->completeOrder(globalOrder, nbrSendOverlap, nbrRecvOverlap, true);
668 if (debug > 1) globalOrder->view("Completed global order");
669 ierr = renumberLocalAdjacencyGraph(mesh, adjGraph, vertexSendOverlap, globalOrder);
670 if (debug > 1) adjGraph->view("Renumbered adjacency graph");
671 }
672 /* Read out adjacency graph */
673 const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::sieve_type> graph = adjBundle->getSieve();
674 const typename Atlas::chart_type& chart = atlas->getChart();
675
676 if (debug > 1) graph->view("Adjacency graph");
677 ierr = PetscMemzero(dnz, numLocalRows/bs * sizeof(PetscInt));CHKERRQ(ierr);
678 ierr = PetscMemzero(onz, numLocalRows/bs * sizeof(PetscInt));CHKERRQ(ierr);
679 for(typename Atlas::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
680 const typename Atlas::point_type& point = *c_iter;
681
682 if (globalOrder->isLocal(point)) {
683 const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::sieve_type::traits::coneSequence>& adj = graph->cone(point);
684 const typename Mesh::order_type::value_type& rIdx = globalOrder->restrictPoint(point)[0];
685 const int row = rIdx.prefix;
686 const int rSize = rIdx.index/bs;
687
688 if ((debug > 1) && ((bs == 1) || rIdx.index%bs)) std::cout << "["<<graph->commRank()<<"]: row "<<row<<": size " << rIdx.index << " bs "<<bs<<std::endl;
689 if (rSize == 0) continue;
690 for(ALE::Mesh<PetscInt,PetscScalar>::sieve_type::traits::coneSequence::iterator v_iter = adj->begin(); v_iter != adj->end(); ++v_iter) {
691 const ALE::Mesh<PetscInt,PetscScalar>::point_type& neighbor = *v_iter;
692 const typename Mesh::order_type::value_type& cIdx = globalOrder->restrictPoint(neighbor)[0];
693 const int& cSize = cIdx.index/bs;
694
695 if ((debug > 1) && ((bs == 1) || cIdx.index%bs)) std::cout << "["<<graph->commRank()<<"]: col "<<cIdx.prefix<<": size " << cIdx.index << " bs "<<bs<<std::endl;
696 if (cSize > 0) {
697 if (globalOrder->isLocal(neighbor)) {
698 for(int r = 0; r < rSize; ++r) {dnz[(row - firstRow)/bs + r] += cSize;}
699 } else {
700 for(int r = 0; r < rSize; ++r) {onz[(row - firstRow)/bs + r] += cSize;}
701 }
702 }
703 }
704 }
705 }
706 if (debug) {
707 for(int r = 0; r < numLocalRows/bs; r++) {
708 std::cout << "["<<rank<<"]: dnz["<<r<<"]: " << dnz[r] << " onz["<<r<<"]: " << onz[r] << std::endl;
709 }
710 }
711 ierr = MatSeqAIJSetPreallocation(A, 0, dnz);CHKERRQ(ierr);
712 ierr = MatMPIAIJSetPreallocation(A, 0, dnz, 0, onz);CHKERRQ(ierr);
713 ierr = MatSeqBAIJSetPreallocation(A, bs, 0, dnz);CHKERRQ(ierr);
714 ierr = MatMPIBAIJSetPreallocation(A, bs, 0, dnz, 0, onz);CHKERRQ(ierr);
715 // TODO: ierr = MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
716 PetscFunctionReturn(0);
717}
718
719#undef __FUNCT__
720#define __FUNCT__ "preallocateOperator"
721template<typename Atlas>
722PetscErrorCode preallocateOperator(const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar> >& mesh, const int bs, const ALE::Obj<Atlas>& atlas, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::order_type>& globalOrder, PetscInt dnz[], PetscInt onz[], Mat A)
723{
724 typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
725 typedef ALE::SieveAlg<FlexMesh> sieve_alg_type;
726 MPI_Comm comm = mesh->comm();
727 const ALE::Obj<FlexMesh> adjBundle = new FlexMesh(comm, mesh->debug());
728 const ALE::Obj<FlexMesh::sieve_type> adjGraph = new FlexMesh::sieve_type(comm, mesh->debug());
729 const bool bigDebug = mesh->debug() > 1;
730 PetscInt numLocalRows, firstRow;
731 ///PetscInt *dnz, *onz;
732 PetscErrorCode ierr;
733
734 PetscFunctionBegin;
735 adjBundle->setSieve(adjGraph);
736 numLocalRows = globalOrder->getLocalSize();
737 firstRow = globalOrder->getGlobalOffsets()[mesh->commRank()];
738 ///ierr = PetscMalloc2(numLocalRows, PetscInt, &dnz, numLocalRows, PetscInt, &onz);CHKERRQ(ierr);
739 /* Create local adjacency graph */
740 /* In general, we need to get FIAT info that attaches dual basis vectors to sieve points */
741 const typename Atlas::chart_type& chart = atlas->getChart();
742
743 for(typename Atlas::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
744 const Obj<typename sieve_alg_type::supportArray>& star = sieve_alg_type::star(mesh, *c_iter);
745
746 for(typename sieve_alg_type::supportArray::const_iterator s_iter = star->begin(); s_iter != star->end(); ++s_iter) {
747 const Obj<typename sieve_alg_type::coneArray>& closure = sieve_alg_type::closure(mesh, *s_iter);
748
749 for(typename sieve_alg_type::coneArray::const_iterator cl_iter = closure->begin(); cl_iter != closure->end(); ++cl_iter) {
750 adjGraph->addCone(*cl_iter, *c_iter);
751 }
752 }
753 }
754 if (bigDebug) adjGraph->view("Adjacency graph");
755 /* Distribute adjacency graph */
756 adjBundle->constructOverlap();
757 typedef typename FlexMesh::sieve_type::point_type point_type;
758 typedef typename FlexMesh::send_overlap_type send_overlap_type;
759 typedef typename FlexMesh::recv_overlap_type recv_overlap_type;
760 typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, point_type> > send_section_type;
761 typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, point_type> > recv_section_type;
762 const Obj<send_overlap_type>& vertexSendOverlap = mesh->getSendOverlap();
763 const Obj<recv_overlap_type>& vertexRecvOverlap = mesh->getRecvOverlap();
764 const Obj<send_overlap_type> nbrSendOverlap = new send_overlap_type(comm, mesh->debug());
765 const Obj<recv_overlap_type> nbrRecvOverlap = new recv_overlap_type(comm, mesh->debug());
766 const Obj<send_section_type> sendSection = new send_section_type(comm, mesh->debug());
767 const Obj<recv_section_type> recvSection = new recv_section_type(comm, sendSection->getTag(), mesh->debug());
768
769 ALE::Distribution<FlexMesh>::coneCompletion(vertexSendOverlap, vertexRecvOverlap, adjBundle, sendSection, recvSection);
770 /* Distribute indices for new points */
771 ///ALE::Distribution<FlexMesh>::updateOverlap(sendSection, recvSection, nbrSendOverlap, nbrRecvOverlap);
772 ALE::Distribution<FlexMesh>::updateOverlap(vertexSendOverlap, vertexRecvOverlap, sendSection, recvSection, nbrSendOverlap, nbrRecvOverlap);
773 mesh->getFactory()->completeOrder(globalOrder, nbrSendOverlap, nbrRecvOverlap, true);
774 /* Read out adjacency graph */
775 const ALE::Obj<FlexMesh::sieve_type> graph = adjBundle->getSieve();
776
777 ierr = PetscMemzero(dnz, numLocalRows/bs * sizeof(PetscInt));CHKERRQ(ierr);
778 ierr = PetscMemzero(onz, numLocalRows/bs * sizeof(PetscInt));CHKERRQ(ierr);
779 for(typename Atlas::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
780 const typename Atlas::point_type& point = *c_iter;
781
782 if (globalOrder->isLocal(point)) {
783 const ALE::Obj<FlexMesh::sieve_type::traits::coneSequence>& adj = graph->cone(point);
784 const FlexMesh::order_type::value_type& rIdx = globalOrder->restrictPoint(point)[0];
785 const int row = rIdx.prefix;
786 const int rSize = rIdx.index/bs;
787
788 if (bigDebug && rIdx.index%bs) std::cout << "["<<graph->commRank()<<"]: row "<<row<<": size " << rIdx.index << " bs "<<bs<<std::endl;
789 if (rSize == 0) continue;
790 for(FlexMesh::sieve_type::traits::coneSequence::iterator v_iter = adj->begin(); v_iter != adj->end(); ++v_iter) {
791 const FlexMesh::point_type& neighbor = *v_iter;
792 const FlexMesh::order_type::value_type& cIdx = globalOrder->restrictPoint(neighbor)[0];
793 const int& cSize = cIdx.index/bs;
794
795 if (bigDebug && cIdx.index%bs) std::cout << "["<<graph->commRank()<<"]: col "<<cIdx.prefix<<": size " << cIdx.index << " bs "<<bs<<std::endl;
796 if (cSize > 0) {
797 if (globalOrder->isLocal(neighbor)) {
798 for(int r = 0; r < rSize; ++r) {dnz[(row - firstRow)/bs + r] += cSize;}
799 } else {
800 for(int r = 0; r < rSize; ++r) {onz[(row - firstRow)/bs + r] += cSize;}
801 }
802 }
803 }
804 }
805 }
806 if (mesh->debug()) {
807 int rank = mesh->commRank();
808 for(int r = 0; r < numLocalRows/bs; r++) {
809 std::cout << "["<<rank<<"]: dnz["<<r<<"]: " << dnz[r] << " onz["<<r<<"]: " << onz[r] << std::endl;
810 }
811 }
812 ierr = MatSeqAIJSetPreallocation(A, 0, dnz);CHKERRQ(ierr);
813 ierr = MatMPIAIJSetPreallocation(A, 0, dnz, 0, onz);CHKERRQ(ierr);
814 ierr = MatSeqBAIJSetPreallocation(A, bs, 0, dnz);CHKERRQ(ierr);
815 ierr = MatMPIBAIJSetPreallocation(A, bs, 0, dnz, 0, onz);CHKERRQ(ierr);
816 ///ierr = PetscFree2(dnz, onz);CHKERRQ(ierr);
817 ///TODO: ierr = MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
818 PetscFunctionReturn(0);
819}
820
821#undef __FUNCT__
822#define __FUNCT__ "preallocateOperator"
823template<typename Atlas>
824PetscErrorCode preallocateOperator(const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar> >& mesh, const int bs, const ALE::Obj<Atlas>& rowAtlas, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::order_type>& rowGlobalOrder, const ALE::Obj<Atlas>& colAtlas, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::order_type>& colGlobalOrder, Mat A)
825{
826 typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
827 typedef ALE::SieveAlg<FlexMesh> sieve_alg_type;
828 MPI_Comm comm = mesh->comm();
829 const ALE::Obj<FlexMesh> adjBundle = new FlexMesh(comm, mesh->debug());
830 const ALE::Obj<FlexMesh::sieve_type> adjGraph = new FlexMesh::sieve_type(comm, mesh->debug());
831 PetscInt numLocalRows, firstRow;
832 PetscInt *dnz, *onz;
833 PetscErrorCode ierr;
834
835 PetscFunctionBegin;
836 adjBundle->setSieve(adjGraph);
837 numLocalRows = rowGlobalOrder->getLocalSize();
838 firstRow = rowGlobalOrder->getGlobalOffsets()[mesh->commRank()];
839 ierr = PetscMalloc2(numLocalRows, PetscInt, &dnz, numLocalRows, PetscInt, &onz);CHKERRQ(ierr);
840 /* Create local adjacency graph */
841 /* In general, we need to get FIAT info that attaches dual basis vectors to sieve points */
842 const typename Atlas::chart_type& chart = rowAtlas->getChart();
843
844 for(typename Atlas::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
845 const Obj<typename sieve_alg_type::supportArray>& star = sieve_alg_type::star(mesh, *c_iter);
846
847 for(typename sieve_alg_type::supportArray::const_iterator s_iter = star->begin(); s_iter != star->end(); ++s_iter) {
848 const Obj<typename sieve_alg_type::coneArray>& closure = sieve_alg_type::closure(mesh, *s_iter);
849
850 for(typename sieve_alg_type::coneArray::const_iterator cl_iter = closure->begin(); cl_iter != closure->end(); ++cl_iter) {
851 adjGraph->addCone(*cl_iter, *c_iter);
852 }
853 }
854 }
855 /* Distribute adjacency graph */
856 adjBundle->constructOverlap();
857 typedef typename FlexMesh::sieve_type::point_type point_type;
858 typedef typename FlexMesh::send_overlap_type send_overlap_type;
859 typedef typename FlexMesh::recv_overlap_type recv_overlap_type;
860 typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, point_type> > send_section_type;
861 typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, point_type> > recv_section_type;
862 const Obj<send_overlap_type>& vertexSendOverlap = mesh->getSendOverlap();
863 const Obj<recv_overlap_type>& vertexRecvOverlap = mesh->getRecvOverlap();
864 const Obj<send_overlap_type> nbrSendOverlap = new send_overlap_type(comm, mesh->debug());
865 const Obj<recv_overlap_type> nbrRecvOverlap = new recv_overlap_type(comm, mesh->debug());
866 const Obj<send_section_type> sendSection = new send_section_type(comm, mesh->debug());
867 const Obj<recv_section_type> recvSection = new recv_section_type(comm, sendSection->getTag(), mesh->debug());
868
869 ALE::Distribution<FlexMesh>::coneCompletion(vertexSendOverlap, vertexRecvOverlap, adjBundle, sendSection, recvSection);
870 /* Distribute indices for new points */
871 ALE::Distribution<FlexMesh>::updateOverlap(vertexSendOverlap, vertexRecvOverlap, sendSection, recvSection, nbrSendOverlap, nbrRecvOverlap);
872 mesh->getFactory()->completeOrder(rowGlobalOrder, nbrSendOverlap, nbrRecvOverlap, true);
873 mesh->getFactory()->completeOrder(colGlobalOrder, nbrSendOverlap, nbrRecvOverlap, true);
874 /* Read out adjacency graph */
875 const ALE::Obj<FlexMesh::sieve_type> graph = adjBundle->getSieve();
876
877 ierr = PetscMemzero(dnz, numLocalRows/bs * sizeof(PetscInt));CHKERRQ(ierr);
878 ierr = PetscMemzero(onz, numLocalRows/bs * sizeof(PetscInt));CHKERRQ(ierr);
879 for(typename Atlas::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
880 const typename Atlas::point_type& point = *c_iter;
881
882 if (rowGlobalOrder->isLocal(point)) {
883 const ALE::Obj<FlexMesh::sieve_type::traits::coneSequence>& adj = graph->cone(point);
884 const FlexMesh::order_type::value_type& rIdx = rowGlobalOrder->restrictPoint(point)[0];
885 const int row = rIdx.prefix;
886 const int rSize = rIdx.index/bs;
887
888 //if (rIdx.index%bs) std::cout << "["<<graph->commRank()<<"]: row "<<row<<": size " << rIdx.index << " bs "<<bs<<std::endl;
889 if (rSize == 0) continue;
890 for(FlexMesh::sieve_type::traits::coneSequence::iterator v_iter = adj->begin(); v_iter != adj->end(); ++v_iter) {
891 const FlexMesh::point_type& neighbor = *v_iter;
892 const FlexMesh::order_type::value_type& cIdx = colGlobalOrder->restrictPoint(neighbor)[0];
893 const int& cSize = cIdx.index/bs;
894
895 //if (cIdx.index%bs) std::cout << "["<<graph->commRank()<<"]: col "<<cIdx.prefix<<": size " << cIdx.index << " bs "<<bs<<std::endl;
896 if (cSize > 0) {
897 if (colGlobalOrder->isLocal(neighbor)) {
898 for(int r = 0; r < rSize; ++r) {dnz[(row - firstRow)/bs + r] += cSize;}
899 } else {
900 for(int r = 0; r < rSize; ++r) {onz[(row - firstRow)/bs + r] += cSize;}
901 }
902 }
903 }
904 }
905 }
906 if (mesh->debug()) {
907 int rank = mesh->commRank();
908 for(int r = 0; r < numLocalRows/bs; r++) {
909 std::cout << "["<<rank<<"]: dnz["<<r<<"]: " << dnz[r] << " onz["<<r<<"]: " << onz[r] << std::endl;
910 }
911 }
912 ierr = MatSeqAIJSetPreallocation(A, 0, dnz);CHKERRQ(ierr);
913 ierr = MatMPIAIJSetPreallocation(A, 0, dnz, 0, onz);CHKERRQ(ierr);
914 ierr = MatSeqBAIJSetPreallocation(A, bs, 0, dnz);CHKERRQ(ierr);
915 ierr = MatMPIBAIJSetPreallocation(A, bs, 0, dnz, 0, onz);CHKERRQ(ierr);
916 ierr = PetscFree2(dnz, onz);CHKERRQ(ierr);
917 ierr = MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
918 PetscFunctionReturn(0);
919}
920
921#undef __FUNCT__
922#define __FUNCT__ "createAllocationVectors"
923template<typename Atlas, typename Order>
924PetscErrorCode createAllocationVectors(const int bs, const ALE::Obj<Atlas>& atlas, const ALE::Obj<Order>& globalOrder, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::sieve_type>& adjGraph, PetscBool isSymmetric, PetscInt dnz[], PetscInt onz[])
925{
926 typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
927 const typename Atlas::chart_type& chart = atlas->getChart();
928 PetscInt numLocalRows = globalOrder->getLocalSize();
929 PetscInt firstRow = globalOrder->getGlobalOffsets()[atlas->commRank()];
930 PetscErrorCode ierr;
931
932 PetscFunctionBegin;
933 ierr = PetscMemzero(dnz, numLocalRows/bs * sizeof(PetscInt));CHKERRQ(ierr);
934 ierr = PetscMemzero(onz, numLocalRows/bs * sizeof(PetscInt));CHKERRQ(ierr);
935 for(typename Atlas::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
936 const typename Atlas::point_type& point = *c_iter;
937
938 if (globalOrder->isLocal(point)) {
939 const ALE::Obj<typename FlexMesh::sieve_type::coneSequence>& adj = adjGraph->cone(point);
940 const typename Order::value_type& rIdx = globalOrder->restrictPoint(point)[0];
941 const int row = rIdx.prefix;
942 const int rSize = rIdx.index/bs;
943
944 if ((atlas->debug() > 1) && ((bs == 1) || (rIdx.index%bs == 0))) std::cout << "["<<adjGraph->commRank()<<"]: row "<<row<<": size " << rIdx.index << " bs "<<bs<<std::endl;
945 if (rSize == 0) continue;
946 for(typename FlexMesh::sieve_type::coneSequence::iterator v_iter = adj->begin(); v_iter != adj->end(); ++v_iter) {
947 const typename Atlas::point_type& neighbor = *v_iter;
948 const typename Order::value_type& cIdx = globalOrder->restrictPoint(neighbor)[0];
949 const int col = cIdx.prefix>=0 ? cIdx.prefix : -(cIdx.prefix+1);
950 const int& cSize = cIdx.index/bs;
951
952 if ((atlas->debug() > 1) && ((bs == 1) || (cIdx.index%bs == 0))) std::cout << "["<<adjGraph->commRank()<<"]: col "<<col<<": size " << cIdx.index << " bs "<<bs<<std::endl;
953 if (cSize > 0) {
954 if (isSymmetric && (col < row)) {
955 if (atlas->debug() > 1) {std::cout << "["<<adjGraph->commRank()<<"]: Rejecting row "<<row<<" col " << col <<std::endl;}
956 continue;
957 }
958 if (globalOrder->isLocal(neighbor)) {
959 for(int r = 0; r < rSize; ++r) {dnz[(row - firstRow)/bs + r] += cSize;}
960 } else {
961 for(int r = 0; r < rSize; ++r) {onz[(row - firstRow)/bs + r] += cSize;}
962 }
963 }
964 }
965 }
966 }
967 PetscFunctionReturn(0);
968}
969
970#undef __FUNCT__
971#define __FUNCT__ "fillMatrixWithZero"
972template<typename Atlas, typename Order>
973PetscErrorCode fillMatrixWithZero(Mat A, const int bs, const ALE::Obj<Atlas>& atlas, const ALE::Obj<Order>& globalOrder, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::sieve_type>& adjGraph, PetscBool isSymmetric, PetscInt dnz[], PetscInt onz[])
974{
975 typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
976 const typename Atlas::chart_type& chart = atlas->getChart();
977 PetscInt numLocalRows = globalOrder->getLocalSize();
978 PetscInt firstRow = globalOrder->getGlobalOffsets()[atlas->commRank()];
979 PetscInt maxRowLen = 0;
980 PetscErrorCode ierr;
981
982 PetscFunctionBegin;
983 for(PetscInt r = 0; r < numLocalRows/bs; ++r) {
984 maxRowLen = std::max(maxRowLen, dnz[r] + onz[r]);
985 }
986 PetscInt *cols = new PetscInt[maxRowLen];
987 PetscScalar *values = new PetscScalar[maxRowLen];
988
989 ierr = PetscMemzero((void *) values, maxRowLen * sizeof(PetscScalar));CHKERRQ(ierr);
990 for(typename Atlas::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
991 const typename Atlas::point_type& point = *c_iter;
992 int rowLen = 0;
993
994 if (globalOrder->isLocal(point)) {
995 const ALE::Obj<typename FlexMesh::sieve_type::coneSequence>& adj = adjGraph->cone(point);
996 const typename Order::value_type& rIdx = globalOrder->restrictPoint(point)[0];
997 const int row = rIdx.prefix;
998 const int rSize = rIdx.index/bs;
999
1000 if (rSize == 0) continue;
1001 for(typename FlexMesh::sieve_type::coneSequence::iterator v_iter = adj->begin(); v_iter != adj->end(); ++v_iter) {
1002 const typename Atlas::point_type& neighbor = *v_iter;
1003 const typename Order::value_type& cIdx = globalOrder->restrictPoint(neighbor)[0];
1004 const int col = cIdx.prefix>=0 ? cIdx.prefix : -(cIdx.prefix+1);
1005 const int& cSize = cIdx.index/bs;
1006
1007 if (cSize > 0) {
1008 if (isSymmetric && (col < row)) {
1009 continue;
1010 }
1011 for(int c = col; c < col+cSize; ++c) {
1012 cols[rowLen++] = c;
1013 }
1014 }
1015 }
1016 for(int r = 0; r < rSize; ++r) {
1017 PetscInt fullRow = row + r;
1018
1019 if (rowLen != dnz[(row - firstRow)/bs+r]+onz[(row - firstRow)/bs+r]) {
1020 SETERRQ5(atlas->comm(), PETSC_ERR_ARG_WRONG, "Invalid row length %d, should be dnz[%d]: %d + onz[%d]: %d", rowLen, (row - firstRow)/bs+r, dnz[(row - firstRow)/bs+r], (row - firstRow)/bs+r, onz[(row - firstRow)/bs+r]);
1021 }
1022 ierr = MatSetValues(A, 1, &fullRow, rowLen, cols, values, INSERT_VALUES);CHKERRQ(ierr);
1023 }
1024 }
1025 }
1026 delete [] cols;
1027 delete [] values;
1028 ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1029 ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1030 PetscFunctionReturn(0);
1031}
1032
1033#undef __FUNCT__
1034#define __FUNCT__ "preallocateOperatorNew"
1035template<typename Mesh, typename Atlas>
1036PetscErrorCode preallocateOperatorNew(const ALE::Obj<Mesh>& mesh, const int bs, const ALE::Obj<Atlas>& atlas, const ALE::Obj<typename Mesh::order_type>& globalOrder, PetscInt dnz[], PetscInt onz[], PetscBool isSymmetric, Mat A, bool fillMatrix = false)
1037{
1038 typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
1039 typedef typename Mesh::sieve_type sieve_type;
1040 typedef typename Mesh::point_type point_type;
1041 typedef typename Mesh::send_overlap_type send_overlap_type;
1042 typedef typename Mesh::recv_overlap_type recv_overlap_type;
1043 const ALE::Obj<FlexMesh::sieve_type> adjGraph = new FlexMesh::sieve_type(mesh->comm(), mesh->debug());
1044 PetscInt numLocalRows = globalOrder->getLocalSize();
1045 PetscInt firstRow = globalOrder->getGlobalOffsets()[mesh->commRank()];
1046 const PetscInt debug = mesh->debug()/3;
1047 PetscErrorCode ierr;
1048
1049 PetscFunctionBegin;
1050 // Create local adjacency graph
1051 if (debug) mesh->view("Input Mesh");
1052 if (debug) globalOrder->view("Initial Global Order");
1053 ierr = createLocalAdjacencyGraph(mesh, atlas, adjGraph);CHKERRQ(ierr);
1054 if (debug) adjGraph->view("Adjacency Graph");
1055 // Complete adjacency graph
1056 typedef ALE::ConeSection<FlexMesh::sieve_type> cones_wrapper_type;
1057 typedef ALE::Section<ALE::Pair<int, point_type>, point_type> cones_type;
1058 Obj<cones_wrapper_type> cones = new cones_wrapper_type(adjGraph);
1059 Obj<cones_type> overlapCones = new cones_type(adjGraph->comm(), adjGraph->debug());
1060 const Obj<send_overlap_type>& sendOverlap = mesh->getSendOverlap();
1061 const Obj<recv_overlap_type>& recvOverlap = mesh->getRecvOverlap();
1062 const Obj<send_overlap_type> nbrSendOverlap = new send_overlap_type(mesh->comm(), mesh->debug());
1063 const Obj<recv_overlap_type> nbrRecvOverlap = new recv_overlap_type(mesh->comm(), mesh->debug());
1064
1065 ALE::Pullback::SimpleCopy::copy(sendOverlap, recvOverlap, cones, overlapCones);
1066 if (debug) overlapCones->view("Overlap Cones");
1067 // Now overlapCones has the neighbors for any point in the overlap, in the remote numbering
1068 // Copy overlaps
1069 {
1070 const Obj<typename send_overlap_type::traits::capSequence> sPoints = sendOverlap->cap();
1071 const typename send_overlap_type::traits::capSequence::iterator sEnd = sPoints->end();
1072
1073 for(typename send_overlap_type::traits::capSequence::iterator p_iter = sPoints->begin(); p_iter != sEnd; ++p_iter) {
1074 const Obj<typename send_overlap_type::supportSequence> support = sendOverlap->support(*p_iter);
1075 const typename send_overlap_type::supportSequence::iterator supEnd = support->end();
1076
1077 for(typename send_overlap_type::supportSequence::iterator s_iter = support->begin(); s_iter != supEnd; ++s_iter) {
1078 nbrSendOverlap->addArrow(*p_iter, *s_iter, s_iter.color());
1079 }
1080 }
1081 const Obj<typename recv_overlap_type::traits::baseSequence> rPoints = recvOverlap->base();
1082 const typename recv_overlap_type::traits::baseSequence::iterator rEnd = rPoints->end();
1083
1084 for(typename recv_overlap_type::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1085 const Obj<typename recv_overlap_type::coneSequence> cone = recvOverlap->cone(*p_iter);
1086 const typename recv_overlap_type::coneSequence::iterator cEnd = cone->end();
1087
1088 for(typename recv_overlap_type::coneSequence::iterator c_iter = cone->begin(); c_iter != cEnd; ++c_iter) {
1089 nbrRecvOverlap->addArrow(*c_iter, *p_iter, c_iter.color());
1090 }
1091 }
1092 }
1093 if (debug) nbrSendOverlap->view("Initial Send Overlap");
1094 if (debug) nbrRecvOverlap->view("Initial Recv Overlap");
1095 // Update neighbor send overlap from local adjacency
1096 typedef typename send_overlap_type::target_type rank_type;
1097 const Obj<typename send_overlap_type::traits::capSequence> sPoints = sendOverlap->cap();
1098 const typename send_overlap_type::traits::capSequence::iterator sEnd = sPoints->end();
1099
1100 for(typename send_overlap_type::traits::capSequence::iterator p_iter = sPoints->begin(); p_iter != sEnd; ++p_iter) {
1101 const point_type& localPoint = *p_iter;
1102 const Obj<typename send_overlap_type::supportSequence>& ranks = sendOverlap->support(localPoint);
1103 const typename send_overlap_type::supportSequence::iterator rEnd = ranks->end();
1104
1105 for(typename send_overlap_type::supportSequence::iterator r_iter = ranks->begin(); r_iter != rEnd; ++r_iter) {
1106 const Obj<typename FlexMesh::sieve_type::coneSequence>& adj = adjGraph->cone(localPoint);
1107 typename FlexMesh::sieve_type::coneSequence::iterator adjEnd = adj->end();
1108
1109 for(typename FlexMesh::sieve_type::coneSequence::iterator c_iter = adj->begin(); c_iter != adjEnd; ++c_iter) {
1110 // Check for interior points
1111 if (!recvOverlap->coneContains(*c_iter, ALE::IsEqual<rank_type>(*r_iter))) {
1112 nbrSendOverlap->addArrow(*c_iter, *r_iter, -1);
1113 }
1114 }
1115 }
1116 }
1117 if (debug) nbrSendOverlap->view("Modified Send Overlap");
1118 // Update neighbor recv overlap and local adjacency
1119 const Obj<typename recv_overlap_type::traits::baseSequence> rPoints = recvOverlap->base();
1120 const typename recv_overlap_type::traits::baseSequence::iterator rEnd = rPoints->end();
1121 point_type maxPoint = std::max(*std::max_element(adjGraph->cap()->begin(), adjGraph->cap()->end()),
1122 *std::max_element(adjGraph->base()->begin(), adjGraph->base()->end())) + 1;
1123
1124 for(typename recv_overlap_type::traits::baseSequence::iterator p_iter = rPoints->begin(); p_iter != rEnd; ++p_iter) {
1125 const point_type& localPoint = *p_iter;
1126 const Obj<typename recv_overlap_type::coneSequence>& ranks = recvOverlap->cone(localPoint);
1127 const typename recv_overlap_type::coneSequence::iterator rEnd = ranks->end();
1128
1129 for(typename recv_overlap_type::coneSequence::iterator r_iter = ranks->begin(); r_iter != rEnd; ++r_iter) {
1130 const int rank = *r_iter;
1131 const point_type& remotePoint = r_iter.color();
1132 const int size = overlapCones->getFiberDimension(typename cones_type::point_type(rank, remotePoint));
1133 const typename cones_type::value_type *values = overlapCones->restrictPoint(typename cones_type::point_type(rank, remotePoint));
1134
1135 for(int i = 0; i < size; ++i) {
1136 // Check for interior point
1137 if (!sendOverlap->cone(rank, values[i])->size()) {
1138 // Check that we have not seen it before
1139 const Obj<typename recv_overlap_type::supportSequence>& newPoints = nbrRecvOverlap->support(rank, values[i]);
1140 point_type newPoint;
1141
1142 if (!newPoints->size()) {
1143 typename Mesh::order_type::value_type value(-1, 0);
1144
1145 newPoint = maxPoint++;
1146 globalOrder->updatePoint(newPoint, &value); // Mark the new point as nonlocal
1147 nbrRecvOverlap->addArrow(rank, newPoint, values[i]);
1148 } else {
1149 newPoint = *newPoints->begin();
1150 }
1151 adjGraph->addArrow(newPoint, localPoint);
1152 adjGraph->addArrow(localPoint, newPoint);
1153 } else {
1154 // Might provide an unknown link for already known point
1155 const point_type oldPoint = *sendOverlap->cone(rank, values[i])->begin();
1156
1157 adjGraph->addArrow(oldPoint, localPoint);
1158 adjGraph->addArrow(localPoint, oldPoint);
1159 }
1160 }
1161 }
1162 }
1163 if (debug) nbrRecvOverlap->view("Modified Recv Overlap");
1164 if (debug) adjGraph->view("Modified Adjacency Graph");
1165 mesh->getFactory()->completeOrder(globalOrder, nbrSendOverlap, nbrRecvOverlap);
1166 if (debug) globalOrder->view("Modified Global Order");
1167 // Read out adjacency graph
1168 const typename Atlas::chart_type& chart = atlas->getChart();
1169
1170 ierr = PetscMemzero(dnz, numLocalRows/bs * sizeof(PetscInt));CHKERRQ(ierr);
1171 ierr = PetscMemzero(onz, numLocalRows/bs * sizeof(PetscInt));CHKERRQ(ierr);
1172 for(typename Atlas::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1173 const typename Atlas::point_type& point = *c_iter;
1174
1175 if (globalOrder->isLocal(point)) {
1176 const ALE::Obj<typename FlexMesh::sieve_type::traits::coneSequence>& adj = adjGraph->cone(point);
1177 const typename Mesh::order_type::value_type& rIdx = globalOrder->restrictPoint(point)[0];
1178 const int row = rIdx.prefix;
1179 const int rSize = rIdx.index/bs;
1180
1181 if ((mesh->debug() > 1) && ((bs == 1) || (rIdx.index%bs == 0))) std::cout << "["<<adjGraph->commRank()<<"]: row "<<row<<": size " << rIdx.index << " bs "<<bs<<std::endl;
1182 if (rSize == 0) continue;
1183 for(typename FlexMesh::sieve_type::traits::coneSequence::iterator v_iter = adj->begin(); v_iter != adj->end(); ++v_iter) {
1184 const typename Mesh::point_type& neighbor = *v_iter;
1185 const typename Mesh::order_type::value_type& cIdx = globalOrder->restrictPoint(neighbor)[0];
1186 const int col = cIdx.prefix>=0 ? cIdx.prefix : -(cIdx.prefix+1);
1187 const int& cSize = cIdx.index/bs;
1188
1189 if ((mesh->debug() > 1) && ((bs == 1) || (cIdx.index%bs == 0))) std::cout << "["<<adjGraph->commRank()<<"]: col "<<col<<": size " << cIdx.index << " bs "<<bs<<std::endl;
1190 if (cSize > 0) {
1191 if (isSymmetric && (col < row)) {
1192 if (mesh->debug() > 1) {std::cout << "["<<adjGraph->commRank()<<"]: Rejecting row "<<row<<" col " << col <<std::endl;}
1193 continue;
1194 }
1195 if (globalOrder->isLocal(neighbor)) {
1196 for(int r = 0; r < rSize; ++r) {dnz[(row - firstRow)/bs + r] += cSize;}
1197 } else {
1198 for(int r = 0; r < rSize; ++r) {onz[(row - firstRow)/bs + r] += cSize;}
1199 }
1200 }
1201 }
1202 }
1203 }
1204 // Set matrix pattern
1205 ierr = MatSeqAIJSetPreallocation(A, 0, dnz);CHKERRQ(ierr);
1206 ierr = MatMPIAIJSetPreallocation(A, 0, dnz, 0, onz);CHKERRQ(ierr);
1207 ierr = MatSeqBAIJSetPreallocation(A, bs, 0, dnz);CHKERRQ(ierr);
1208 ierr = MatMPIBAIJSetPreallocation(A, bs, 0, dnz, 0, onz);CHKERRQ(ierr);
1209 ierr = MatSeqSBAIJSetPreallocation(A, bs, 0, dnz);CHKERRQ(ierr);
1210 ierr = MatMPISBAIJSetPreallocation(A, bs, 0, dnz, 0, onz);CHKERRQ(ierr);
1211 ierr = MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
1212 // Fill matrix with zeros
1213 if (fillMatrix) {
1214 int maxRowLen = 0;
1215 for(int r = 0; r < numLocalRows/bs; ++r) {
1216 maxRowLen = std::max(maxRowLen, dnz[r] + onz[r]);
1217 }
1218 PetscInt *cols = new PetscInt[maxRowLen];
1219 PetscScalar *values = new PetscScalar[maxRowLen];
1220
1221 ierr = PetscMemzero((void *) values, maxRowLen * sizeof(PetscScalar));CHKERRQ(ierr);
1222 for(typename Atlas::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
1223 const typename Atlas::point_type& point = *c_iter;
1224 int rowLen = 0;
1225
1226 if (globalOrder->isLocal(point)) {
1227 const ALE::Obj<typename FlexMesh::sieve_type::traits::coneSequence>& adj = adjGraph->cone(point);
1228 const typename Mesh::order_type::value_type& rIdx = globalOrder->restrictPoint(point)[0];
1229 const int row = rIdx.prefix;
1230 const int rSize = rIdx.index/bs;
1231
1232 if (rSize == 0) continue;
1233 for(typename FlexMesh::sieve_type::traits::coneSequence::iterator v_iter = adj->begin(); v_iter != adj->end(); ++v_iter) {
1234 const typename Mesh::point_type& neighbor = *v_iter;
1235 const typename Mesh::order_type::value_type& cIdx = globalOrder->restrictPoint(neighbor)[0];
1236 const int col = cIdx.prefix>=0 ? cIdx.prefix : -(cIdx.prefix+1);
1237 const int& cSize = cIdx.index/bs;
1238
1239 if (cSize > 0) {
1240 if (isSymmetric && (col < row)) {
1241 continue;
1242 }
1243 for(int c = col; c < col+cSize; ++c) {
1244 cols[rowLen++] = c;
1245 }
1246 }
1247 }
1248 for(int r = 0; r < rSize; ++r) {
1249 PetscInt fullRow = row + r;
1250
1251 if (rowLen != dnz[(row - firstRow)/bs+r]+onz[(row - firstRow)/bs+r]) {
1252 SETERRQ5(mesh->comm(), PETSC_ERR_ARG_WRONG, "Invalid row length %d, should be dnz[%d]: %d + onz[%d]: %d", rowLen, (row - firstRow)/bs+r, dnz[(row - firstRow)/bs+r], (row - firstRow)/bs+r, onz[(row - firstRow)/bs+r]);
1253 }
1254 ierr = MatSetValues(A, 1, &fullRow, rowLen, cols, values, INSERT_VALUES);CHKERRQ(ierr);
1255 }
1256 }
1257 }
1258 delete [] cols;
1259 delete [] values;
1260 ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1261 ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
1262 }
1263 PetscFunctionReturn(0);
1264}
1265
1266#undef __FUNCT__
1267#define __FUNCT__ "preallocateOperatorNewOverlap"
1268 template<typename Mesh, typename Atlas>
1269PetscErrorCode preallocateOperatorNewOverlap(const ALE::Obj<Mesh>& mesh, const int bs, const ALE::Obj<Atlas>& atlas, const ALE::Obj<typename Mesh::order_type>& globalOrder, PetscInt dnz[], PetscInt onz[], PetscBool isSymmetric, Mat A, bool fillMatrix = false)
1270{
1271 typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
1272 typedef typename Mesh::sieve_type sieve_type;
1273 typedef typename Mesh::point_type point_type;
1274 typedef typename Mesh::send_overlap_type send_overlap_type;
1275 typedef typename Mesh::recv_overlap_type recv_overlap_type;
1276 const ALE::Obj<FlexMesh::sieve_type> adjGraph = new FlexMesh::sieve_type(mesh->comm(), mesh->debug());
1277 const PetscInt debug = mesh->debug()/3;
1278 PetscErrorCode ierr;
1279
1280 PetscFunctionBegin;
1281 // Create local adjacency graph
1282 if (debug) mesh->view("Input Mesh");
1283 if (debug) globalOrder->view("Initial Global Order");
1284 ierr = createLocalAdjacencyGraph(mesh, atlas, adjGraph);CHKERRQ(ierr);
1285 if (debug) adjGraph->view("Adjacency Graph");
1286 // Complete adjacency graph
1287 typedef ALE::ConeSection<FlexMesh::sieve_type> cones_wrapper_type;
1288 typedef ALE::Section<ALE::Pair<int, point_type>, point_type> cones_type;
1289 Obj<cones_wrapper_type> cones = new cones_wrapper_type(adjGraph);
1290 Obj<cones_type> overlapCones = new cones_type(adjGraph->comm(), adjGraph->debug());
1291 const Obj<send_overlap_type>& sendOverlap = mesh->getSendOverlap();
1292 const Obj<recv_overlap_type>& recvOverlap = mesh->getRecvOverlap();
1293 const Obj<send_overlap_type> nbrSendOverlap = new send_overlap_type(mesh->comm(), mesh->debug());
1294 const Obj<recv_overlap_type> nbrRecvOverlap = new recv_overlap_type(mesh->comm(), mesh->debug());
1295
1296 // Now overlapCones will have the neighbors for any point in the overlap, in the remote numbering
1297 ALE::Pullback::SimpleCopy::copy(sendOverlap, recvOverlap, cones, overlapCones);
1298 if (debug) overlapCones->view("Overlap Cones");
1299 // TODO Copy overlaps
1300 sendOverlap->copy(nbrSendOverlap.ptr());
1301 recvOverlap->copy(nbrRecvOverlap.ptr());
1302 if (debug) nbrSendOverlap->view("Initial Send Overlap");
1303 if (debug) nbrRecvOverlap->view("Initial Recv Overlap");
1304 // TODO Update neighbor send overlap from local adjacency
1305 // For each localPoint in sendOverlap
1306 // For each rank receiving this point
1307 // For each adjPoint in adjGraph->cone(point)
1308 // If recvOverlap does not contain an arrow (rank, adjPoint, *), meaning the point is not interior to the domain
1309 // nbrSendOverlap->addArrow(adjPoint, rank, -1)
1310 const typename send_overlap_type::baseSequence::iterator sBegin = sendOverlap->baseBegin();
1311 const typename send_overlap_type::baseSequence::iterator sEnd = sendOverlap->baseEnd();
1312
1313 for(typename send_overlap_type::baseSequence::iterator r_iter = sBegin; r_iter != sEnd; ++r_iter) {
1314 const typename send_overlap_type::target_type rank = *r_iter;
1315 const typename send_overlap_type::coneSequence::iterator pBegin = sendOverlap->coneBegin(*r_iter);
1316 const typename send_overlap_type::coneSequence::iterator pEnd = sendOverlap->coneEnd(*r_iter);
1317
1318 for(typename send_overlap_type::coneSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
1319 const typename send_overlap_type::source_type localPoint = *p_iter;
1320 const typename FlexMesh::sieve_type::coneSequence::iterator adjBegin = adjGraph->cone(localPoint)->begin();
1321 const typename FlexMesh::sieve_type::coneSequence::iterator adjEnd = adjGraph->cone(localPoint)->end();
1322
1323 for(typename FlexMesh::sieve_type::coneSequence::iterator a_iter = adjBegin; a_iter != adjEnd; ++a_iter) {
1324 const typename FlexMesh::sieve_type::coneSequence::iterator::value_type adjPoint = *a_iter;
1325
1326 // Deal with duplication at the assembly stage
1327 nbrSendOverlap->addArrow(adjPoint, rank, -1);
1328 }
1329 }
1330 }
1331 nbrSendOverlap->assemble();
1332 nbrSendOverlap->assemblePoints();
1333 if (debug) nbrSendOverlap->view("Modified Send Overlap");
1334 // Let maxPoint be the first point not contained in adjGraph
1335 point_type maxPoint = std::max(*std::max_element(adjGraph->cap()->begin(), adjGraph->cap()->end()),
1336 *std::max_element(adjGraph->base()->begin(), adjGraph->base()->end())) + 1;
1337 // TODO Update neighbor recv overlap and local adjacency
1338 // For each point in recvOverlap
1339 // For each rank sending this point
1340 // For each adjPoint in the overlap cone from adjGraph for this point
1341 // If adjPoint is interior, meaning sendOverlap has no arrow (rank, *, adjPoint) CAN THIS EVER HAPPEN???
1342 // If nbrRevOverlap has arrow (rank, newPoint, adjPoint)
1343 // Let newPoint = maxPoint, increment maxPoint
1344 // Add arrows (point, newPoint) and (newPoint, point) to adjGraph
1345 // Else
1346 // Add arrows (point, newPoint) and (newPoint, point) to adjGraph
1347 // Else
1348 // Why would we see a new connection for an old point??? Need an example
1349 // We have the arrow (rank, oldPoint, adjPoint)
1350 // Add arrows (point, oldPoint) and (oldPoint, point) to adjGraph
1351 const typename recv_overlap_type::capSequence::iterator rBegin = recvOverlap->capBegin();
1352 const typename recv_overlap_type::capSequence::iterator rEnd = recvOverlap->capEnd();
1353
1354 for(typename recv_overlap_type::capSequence::iterator r_iter = rBegin; r_iter != rEnd; ++r_iter) {
1355 const int rank = *r_iter;
1356 const typename recv_overlap_type::supportSequence::iterator pBegin = recvOverlap->supportBegin(*r_iter);
1357 const typename recv_overlap_type::supportSequence::iterator pEnd = recvOverlap->supportEnd(*r_iter);
1358
1359 for(typename recv_overlap_type::supportSequence::iterator p_iter = pBegin; p_iter != pEnd; ++p_iter) {
1360 const point_type& localPoint = *p_iter;
1361 const point_type& remotePoint = p_iter.color();
1362 const int size = overlapCones->getFiberDimension(typename cones_type::point_type(rank, remotePoint));
1363 const typename cones_type::value_type *values = overlapCones->restrictPoint(typename cones_type::point_type(rank, remotePoint));
1364
1365 for(int i = 0; i < size; ++i) {
1366 const typename recv_overlap_type::supportSequence::iterator newPointsBegin = nbrRecvOverlap->supportBegin(rank, values[i]);
1367 const int numNewPoints = nbrRecvOverlap->getSupportSize(rank, values[i]);
1368 point_type newPoint;
1369
1370 if (!numNewPoints) {
1371 typename Mesh::order_type::value_type value(-1, 0);
1372
1373 newPoint = maxPoint++;
1374 globalOrder->updatePoint(newPoint, &value); // Mark the new point as nonlocal
1375 nbrRecvOverlap->addArrow(rank, newPoint, values[i]);
1376 } else {
1377 newPoint = *newPointsBegin;
1378 }
1379 adjGraph->addArrow(newPoint, localPoint);
1380 adjGraph->addArrow(localPoint, newPoint);
1381 }
1382 }
1383 }
1384 nbrRecvOverlap->assemble();
1385 nbrRecvOverlap->assemblePoints();
1386 if (debug) nbrRecvOverlap->view("Modified Recv Overlap");
1387 if (debug) adjGraph->view("Modified Adjacency Graph");
1388 // Update global order
1389 mesh->getFactory()->completeOrder(globalOrder, nbrSendOverlap, nbrRecvOverlap);
1390 if (debug) globalOrder->view("Modified Global Order");
1391 // Read out adjacency graph
1392 ierr = createAllocationVectors(bs, atlas, globalOrder, adjGraph, isSymmetric, dnz, onz);
1393 // Set matrix pattern
1394 ierr = MatSeqAIJSetPreallocation(A, 0, dnz);CHKERRQ(ierr);
1395 ierr = MatMPIAIJSetPreallocation(A, 0, dnz, 0, onz);CHKERRQ(ierr);
1396 ierr = MatSeqBAIJSetPreallocation(A, bs, 0, dnz);CHKERRQ(ierr);
1397 ierr = MatMPIBAIJSetPreallocation(A, bs, 0, dnz, 0, onz);CHKERRQ(ierr);
1398 ierr = MatSeqSBAIJSetPreallocation(A, bs, 0, dnz);CHKERRQ(ierr);
1399 ierr = MatMPISBAIJSetPreallocation(A, bs, 0, dnz, 0, onz);CHKERRQ(ierr);
1400 ierr = MatSetUp(A);
1401 ierr = MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
1402 // Fill matrix with zeros
1403 if (fillMatrix) {ierr = fillMatrixWithZero(A, bs, atlas, globalOrder, adjGraph, isSymmetric, dnz, onz);CHKERRQ(ierr);}
1404 PetscFunctionReturn(0);
1405}
1406
1407template<typename Mesh, typename Atlas>
1408PetscErrorCode preallocateOperatorI(const ALE::Obj<Mesh>& mesh, const int bs, const ALE::Obj<Atlas>& atlas, const ALE::Obj<typename Mesh::order_type>& globalOrder, PetscInt dnz[], PetscInt onz[], PetscBool isSymmetric, Mat A)
1409{
1410 typedef typename Mesh::sieve_type sieve_type;
1411 typedef typename Mesh::point_type point_type;
1412 typedef typename Mesh::send_overlap_type send_overlap_type;
1413 typedef typename Mesh::recv_overlap_type recv_overlap_type;
1414 const ALE::Obj<typename Mesh::sieve_type> adjGraph = new typename Mesh::sieve_type(mesh->comm(), mesh->debug());
1415 PetscInt numLocalRows = globalOrder->getLocalSize();
1416 PetscInt firstRow = globalOrder->getGlobalOffsets()[mesh->commRank()];
1417 const PetscInt debug = 0;
1418 PetscErrorCode ierr;
1419
1420 PetscFunctionBegin;
1421 // Create local adjacency graph
1422 if (debug) mesh->view("Input Mesh");
1423 if (debug) globalOrder->view("Initial Global Order");
1424 adjGraph->setChart(mesh);
1425 ierr = createLocalAdjacencyGraphI(mesh, atlas, adjGraph);CHKERRQ(ierr);
1426 if (debug) adjGraph->view("Adjacency Graph");
1427
1428 // Will have to reallocate() adjGraph after adding arrows
1429
1430 // Rewrite read out from adjGraph to use visitors
1431
1432 // Set matrix pattern
1433 ierr = MatSeqAIJSetPreallocation(A, 0, dnz);CHKERRQ(ierr);
1434 ierr = MatMPIAIJSetPreallocation(A, 0, dnz, 0, onz);CHKERRQ(ierr);
1435 ierr = MatSeqBAIJSetPreallocation(A, bs, 0, dnz);CHKERRQ(ierr);
1436 ierr = MatMPIBAIJSetPreallocation(A, bs, 0, dnz, 0, onz);CHKERRQ(ierr);
1437 ierr = MatSeqSBAIJSetPreallocation(A, bs, 0, dnz);CHKERRQ(ierr);
1438 ierr = MatMPISBAIJSetPreallocation(A, bs, 0, dnz, 0, onz);CHKERRQ(ierr);
1439 ierr = MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
1440 PetscFunctionReturn(0);
1441}
1442
1443#undef __FUNCT__
1444#define __FUNCT__ "updateOperator"
1445template<typename Sieve, typename Visitor>
1446PetscErrorCode updateOperator(Mat A, const Sieve& sieve, Visitor& iV, const PETSC_MESH_TYPE::point_type& e, PetscScalar array[], InsertMode mode)
1447{
1448 PetscFunctionBegin;
1449 ALE::ISieveTraversal<Sieve>::orientedClosure(sieve, e, iV);
1450 const PetscInt *indices = iV.getValues();
1451 const int numIndices = iV.getSize();
1452 PetscErrorCode ierr;
1453
1454 ierr = PetscLogEventBegin(DMMesh_updateOperator,0,0,0,0);CHKERRQ(ierr);
1455 if (sieve.debug()) {
1456 ierr = PetscPrintf(PETSC_COMM_SELF, "[%d]mat for element %d\n", sieve.commRank(), e);CHKERRQ(ierr);
1457 for(int i = 0; i < numIndices; i++) {
1458 ierr = PetscPrintf(PETSC_COMM_SELF, "[%d]mat indices[%d] = %d\n", sieve.commRank(), i, indices[i]);CHKERRQ(ierr);
1459 }
1460 for(int i = 0; i < numIndices; i++) {
1461 ierr = PetscPrintf(PETSC_COMM_SELF, "[%d]", sieve.commRank());CHKERRQ(ierr);
1462 for(int j = 0; j < numIndices; j++) {
1463#ifdef PETSC_USE_COMPLEX
1464 ierr = PetscPrintf(PETSC_COMM_SELF, " (%g,%g)", PetscRealPart(array[i*numIndices+j]), PetscImaginaryPart(array[i*numIndices+j]));CHKERRQ(ierr);
1465#else
1466 ierr = PetscPrintf(PETSC_COMM_SELF, " %g", array[i*numIndices+j]);CHKERRQ(ierr);
1467#endif
1468 }
1469 ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr);
1470 }
1471 }
1472 ierr = MatSetValues(A, numIndices, indices, numIndices, indices, array, mode);
1473 if (ierr) {
1474 PetscErrorCode ierr2;
1475 ierr2 = PetscPrintf(PETSC_COMM_SELF, "[%d]ERROR in updateOperator: point %d\n", sieve.commRank(), e);CHKERRQ(ierr2);
1476 for(int i = 0; i < numIndices; i++) {
1477 ierr2 = PetscPrintf(PETSC_COMM_SELF, "[%d]mat indices[%d] = %d\n", sieve.commRank(), i, indices[i]);CHKERRQ(ierr2);
1478 }
1479 CHKERRQ(ierr);
1480 }
1481 ierr = PetscLogEventEnd(DMMesh_updateOperator,0,0,0,0);CHKERRQ(ierr);
1482 PetscFunctionReturn(0);
1483}
1484
1485#undef __FUNCT__
1486#define __FUNCT__ "updateOperator"
1487template<typename Sieve, typename Visitor>
1488PetscErrorCode updateOperator(Mat A, const Sieve& rowSieve, Visitor& iVr, const PETSC_MESH_TYPE::point_type& rowE, const Sieve& colSieve, Visitor& iVc, const PETSC_MESH_TYPE::point_type& colE, PetscScalar array[], InsertMode mode)
1489{
1490 PetscFunctionBegin;
1491 ALE::ISieveTraversal<Sieve>::orientedClosure(rowSieve, rowE, iVr);
1492 ALE::ISieveTraversal<Sieve>::orientedClosure(colSieve, colE, iVc);
1493 const PetscInt *rowIndices = iVr.getValues();
1494 const int numRowIndices = iVr.getSize();
1495 const PetscInt *colIndices = iVc.getValues();
1496 const int numColIndices = iVc.getSize();
1497 PetscErrorCode ierr;
1498
1499 ierr = PetscLogEventBegin(DMMesh_updateOperator,0,0,0,0);CHKERRQ(ierr);
1500 if (rowSieve.debug()) {
1501 ierr = PetscPrintf(PETSC_COMM_SELF, "[%d]mat for element %d,%d\n", rowSieve.commRank(), rowE, colE);CHKERRQ(ierr);
1502 for(int i = 0; i < numRowIndices; i++) {
1503 ierr = PetscPrintf(PETSC_COMM_SELF, "[%d]mat row indices[%d] = %d\n", rowSieve.commRank(), i, rowIndices[i]);CHKERRQ(ierr);
1504 }
1505 for(int i = 0; i < numColIndices; i++) {
1506 ierr = PetscPrintf(PETSC_COMM_SELF, "[%d]mat col indices[%d] = %d\n", rowSieve.commRank(), i, colIndices[i]);CHKERRQ(ierr);
1507 }
1508 for(int i = 0; i < numRowIndices; i++) {
1509 ierr = PetscPrintf(PETSC_COMM_SELF, "[%d]", rowSieve.commRank());CHKERRQ(ierr);
1510 for(int j = 0; j < numColIndices; j++) {
1511#ifdef PETSC_USE_COMPLEX
1512 ierr = PetscPrintf(PETSC_COMM_SELF, " (%g,%g)", PetscRealPart(array[i*numColIndices+j]), PetscImaginaryPart(array[i*numColIndices+j]));CHKERRQ(ierr);
1513#else
1514 ierr = PetscPrintf(PETSC_COMM_SELF, " %g", array[i*numColIndices+j]);CHKERRQ(ierr);
1515#endif
1516 }
1517 ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr);
1518 }
1519 }
1520 ierr = MatSetValues(A, numRowIndices, rowIndices, numColIndices, colIndices, array, mode);
1521 if (ierr) {
1522 PetscErrorCode ierr2;
1523 ierr2 = PetscPrintf(PETSC_COMM_SELF, "[%d]ERROR in updateOperator: point %d,%d\n", rowSieve.commRank(), rowE, colE);CHKERRQ(ierr2);
1524 for(int i = 0; i < numRowIndices; i++) {
1525 ierr2 = PetscPrintf(PETSC_COMM_SELF, "[%d]mat row indices[%d] = %d\n", rowSieve.commRank(), i, rowIndices[i]);CHKERRQ(ierr2);
1526 }
1527 for(int i = 0; i < numColIndices; i++) {
1528 ierr2 = PetscPrintf(PETSC_COMM_SELF, "[%d]mat col indices[%d] = %d\n", rowSieve.commRank(), i, colIndices[i]);CHKERRQ(ierr2);
1529 }
1530 CHKERRQ(ierr);
1531 }
1532 ierr = PetscLogEventEnd(DMMesh_updateOperator,0,0,0,0);CHKERRQ(ierr);
1533 PetscFunctionReturn(0);
1534}
1535
1536#endif // __PETSCDMMESH_HH
Note: See TracBrowser for help on using the repository browser.