/[escript]/trunk/ripley/src/Rectangle.cpp
ViewVC logotype

Contents of /trunk/ripley/src/Rectangle.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3699 - (show annotations)
Thu Dec 1 22:59:42 2011 UTC (7 years, 11 months ago) by caltinay
Original Path: branches/ripleygmg_from_3668/ripley/src/Rectangle.cpp
File size: 29713 byte(s)
getPattern() for Rectangle. Couple node numbering is likely to change once
we get to the PDE assembly.

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2011 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14 #include <ripley/Rectangle.h>
15 extern "C" {
16 #include "paso/SystemMatrixPattern.h"
17 }
18
19 #if USE_SILO
20 #include <silo.h>
21 #ifdef ESYS_MPI
22 #include <pmpio.h>
23 #endif
24 #endif
25
26 #include <iomanip>
27
28 using namespace std;
29
30 namespace ripley {
31
32 Rectangle::Rectangle(int n0, int n1, double l0, double l1, int d0, int d1) :
33 RipleyDomain(2),
34 m_gNE0(n0),
35 m_gNE1(n1),
36 m_l0(l0),
37 m_l1(l1),
38 m_NX(d0),
39 m_NY(d1)
40 {
41 // ensure number of subdivisions is valid and nodes can be distributed
42 // among number of ranks
43 if (m_NX*m_NY != m_mpiInfo->size)
44 throw RipleyException("Invalid number of spatial subdivisions");
45
46 if (n0%m_NX > 0 || n1%m_NY > 0)
47 throw RipleyException("Number of elements must be separable into number of ranks in each dimension");
48
49 // local number of elements
50 m_NE0 = n0/m_NX;
51 m_NE1 = n1/m_NY;
52 // local number of nodes (not necessarily owned)
53 m_N0 = m_NE0+1;
54 m_N1 = m_NE1+1;
55 // bottom-left node is at (offset0,offset1) in global mesh
56 m_offset0 = m_NE0*(m_mpiInfo->rank%m_NX);
57 m_offset1 = m_NE1*(m_mpiInfo->rank/m_NX);
58 populateSampleIds();
59 }
60
61 Rectangle::~Rectangle()
62 {
63 }
64
65 string Rectangle::getDescription() const
66 {
67 return "ripley::Rectangle";
68 }
69
70 bool Rectangle::operator==(const AbstractDomain& other) const
71 {
72 if (dynamic_cast<const Rectangle*>(&other))
73 return this==&other;
74
75 return false;
76 }
77
78 void Rectangle::dump(const string& fileName) const
79 {
80 #if USE_SILO
81 string fn(fileName);
82 if (fileName.length() < 6 || fileName.compare(fileName.length()-5, 5, ".silo") != 0) {
83 fn+=".silo";
84 }
85
86 const int NUM_SILO_FILES = 1;
87 const char* blockDirFmt = "/block%04d";
88 int driver=DB_HDF5;
89 string siloPath;
90 DBfile* dbfile = NULL;
91
92 #ifdef ESYS_MPI
93 PMPIO_baton_t* baton = NULL;
94 #endif
95
96 if (m_mpiInfo->size > 1) {
97 #ifdef ESYS_MPI
98 baton = PMPIO_Init(NUM_SILO_FILES, PMPIO_WRITE, m_mpiInfo->comm,
99 0x1337, PMPIO_DefaultCreate, PMPIO_DefaultOpen,
100 PMPIO_DefaultClose, (void*)&driver);
101 // try the fallback driver in case of error
102 if (!baton && driver != DB_PDB) {
103 driver = DB_PDB;
104 baton = PMPIO_Init(NUM_SILO_FILES, PMPIO_WRITE, m_mpiInfo->comm,
105 0x1338, PMPIO_DefaultCreate, PMPIO_DefaultOpen,
106 PMPIO_DefaultClose, (void*)&driver);
107 }
108 if (baton) {
109 char str[64];
110 snprintf(str, 64, blockDirFmt, PMPIO_RankInGroup(baton, m_mpiInfo->rank));
111 siloPath = str;
112 dbfile = (DBfile*) PMPIO_WaitForBaton(baton, fn.c_str(), siloPath.c_str());
113 }
114 #endif
115 } else {
116 dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,
117 getDescription().c_str(), driver);
118 // try the fallback driver in case of error
119 if (!dbfile && driver != DB_PDB) {
120 driver = DB_PDB;
121 dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,
122 getDescription().c_str(), driver);
123 }
124 }
125
126 if (!dbfile)
127 throw RipleyException("dump: Could not create Silo file");
128
129 /*
130 if (driver==DB_HDF5) {
131 // gzip level 1 already provides good compression with minimal
132 // performance penalty. Some tests showed that gzip levels >3 performed
133 // rather badly on escript data both in terms of time and space
134 DBSetCompression("ERRMODE=FALLBACK METHOD=GZIP LEVEL=1");
135 }
136 */
137
138 boost::scoped_ptr<double> x(new double[m_N0]);
139 boost::scoped_ptr<double> y(new double[m_N1]);
140 double* coords[2] = { x.get(), y.get() };
141 pair<double,double> xdx = getFirstCoordAndSpacing(0);
142 pair<double,double> ydy = getFirstCoordAndSpacing(1);
143 #pragma omp parallel
144 {
145 #pragma omp for
146 for (dim_t i0 = 0; i0 < m_N0; i0++) {
147 coords[0][i0]=xdx.first+i0*xdx.second;
148 }
149 #pragma omp for
150 for (dim_t i1 = 0; i1 < m_N1; i1++) {
151 coords[1][i1]=ydy.first+i1*ydy.second;
152 }
153 }
154 IndexVector dims = getNumNodesPerDim();
155
156 // write mesh
157 DBPutQuadmesh(dbfile, "mesh", NULL, coords, &dims[0], 2, DB_DOUBLE,
158 DB_COLLINEAR, NULL);
159
160 // write node ids
161 DBPutQuadvar1(dbfile, "nodeId", "mesh", (void*)&m_nodeId[0], &dims[0], 2,
162 NULL, 0, DB_INT, DB_NODECENT, NULL);
163
164 // write element ids
165 dims = getNumElementsPerDim();
166 DBPutQuadvar1(dbfile, "elementId", "mesh", (void*)&m_elementId[0],
167 &dims[0], 2, NULL, 0, DB_INT, DB_ZONECENT, NULL);
168
169 // rank 0 writes multimesh and multivar
170 if (m_mpiInfo->rank == 0) {
171 vector<string> tempstrings;
172 vector<char*> names;
173 for (dim_t i=0; i<m_mpiInfo->size; i++) {
174 stringstream path;
175 path << "/block" << setw(4) << setfill('0') << right << i << "/mesh";
176 tempstrings.push_back(path.str());
177 names.push_back((char*)tempstrings.back().c_str());
178 }
179 vector<int> types(m_mpiInfo->size, DB_QUAD_RECT);
180 DBSetDir(dbfile, "/");
181 DBPutMultimesh(dbfile, "multimesh", m_mpiInfo->size, &names[0],
182 &types[0], NULL);
183 tempstrings.clear();
184 names.clear();
185 for (dim_t i=0; i<m_mpiInfo->size; i++) {
186 stringstream path;
187 path << "/block" << setw(4) << setfill('0') << right << i << "/nodeId";
188 tempstrings.push_back(path.str());
189 names.push_back((char*)tempstrings.back().c_str());
190 }
191 types.assign(m_mpiInfo->size, DB_QUADVAR);
192 DBPutMultivar(dbfile, "nodeId", m_mpiInfo->size, &names[0],
193 &types[0], NULL);
194 tempstrings.clear();
195 names.clear();
196 for (dim_t i=0; i<m_mpiInfo->size; i++) {
197 stringstream path;
198 path << "/block" << setw(4) << setfill('0') << right << i << "/elementId";
199 tempstrings.push_back(path.str());
200 names.push_back((char*)tempstrings.back().c_str());
201 }
202 DBPutMultivar(dbfile, "elementId", m_mpiInfo->size, &names[0],
203 &types[0], NULL);
204 }
205
206 if (m_mpiInfo->size > 1) {
207 #ifdef ESYS_MPI
208 PMPIO_HandOffBaton(baton, dbfile);
209 PMPIO_Finish(baton);
210 #endif
211 } else {
212 DBClose(dbfile);
213 }
214
215 #else // USE_SILO
216 throw RipleyException("dump(): no Silo support");
217 #endif
218 }
219
220 const int* Rectangle::borrowSampleReferenceIDs(int fsType) const
221 {
222 switch (fsType) {
223 case Nodes:
224 return &m_nodeId[0];
225 case Elements:
226 return &m_elementId[0];
227 case FaceElements:
228 return &m_faceId[0];
229 default:
230 break;
231 }
232
233 stringstream msg;
234 msg << "borrowSampleReferenceIDs() not implemented for function space type "
235 << fsType;
236 throw RipleyException(msg.str());
237 }
238
239 bool Rectangle::ownSample(int fsCode, index_t id) const
240 {
241 #ifdef ESYS_MPI
242 if (fsCode == Nodes) {
243 const index_t myFirst=m_nodeDistribution[m_mpiInfo->rank];
244 const index_t myLast=m_nodeDistribution[m_mpiInfo->rank+1]-1;
245 return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);
246 } else
247 throw RipleyException("ownSample() only implemented for Nodes");
248 #else
249 return true;
250 #endif
251 }
252
253 void Rectangle::interpolateOnDomain(escript::Data& target,
254 const escript::Data& in) const
255 {
256 const Rectangle& inDomain=dynamic_cast<const Rectangle&>(*(in.getFunctionSpace().getDomain()));
257 const Rectangle& targetDomain=dynamic_cast<const Rectangle&>(*(target.getFunctionSpace().getDomain()));
258 if (inDomain != *this)
259 throw RipleyException("Illegal domain of interpolant");
260 if (targetDomain != *this)
261 throw RipleyException("Illegal domain of interpolation target");
262
263 stringstream msg;
264 msg << "interpolateOnDomain() not implemented for function space "
265 << in.getFunctionSpace().getTypeCode() << " -> "
266 << target.getFunctionSpace().getTypeCode();
267
268 switch (in.getFunctionSpace().getTypeCode()) {
269 case Nodes:
270 case DegreesOfFreedom:
271 switch (target.getFunctionSpace().getTypeCode()) {
272 case DegreesOfFreedom:
273 target=in;
274 break;
275
276 case Elements:
277 {
278 const double tmp0_2 = 0.62200846792814621559;
279 const double tmp0_1 = 0.044658198738520451079;
280 const double tmp0_0 = 0.16666666666666666667;
281 const dim_t numComp = in.getDataPointSize();
282 escript::Data* inn=const_cast<escript::Data*>(&in);
283 #pragma omp parallel for
284 for (index_t k1=0; k1 < m_NE1; ++k1) {
285 for (index_t k0=0; k0 < m_NE0; ++k0) {
286 const register double* f_10 = inn->getSampleDataRO(INDEX2(k0+1,k1, m_N0));
287 const register double* f_11 = inn->getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
288 const register double* f_01 = inn->getSampleDataRO(INDEX2(k0,k1+1, m_N0));
289 const register double* f_00 = inn->getSampleDataRO(INDEX2(k0,k1, m_N0));
290 double* o = target.getSampleDataRW(INDEX2(k0,k1,m_NE0));
291 for (index_t i=0; i < numComp; ++i) {
292 o[INDEX2(i,numComp,0)] = f_00[i]*tmp0_2 + f_11[i]*tmp0_1 + tmp0_0*(f_01[i] + f_10[i]);
293 o[INDEX2(i,numComp,1)] = f_01[i]*tmp0_1 + f_10[i]*tmp0_2 + tmp0_0*(f_00[i] + f_11[i]);
294 o[INDEX2(i,numComp,2)] = f_01[i]*tmp0_2 + f_10[i]*tmp0_1 + tmp0_0*(f_00[i] + f_11[i]);
295 o[INDEX2(i,numComp,3)] = f_00[i]*tmp0_1 + f_11[i]*tmp0_2 + tmp0_0*(f_01[i] + f_10[i]);
296 } // close component loop i
297 } // close k0 loop
298 } // close k1 loop
299 break;
300 }
301
302 default:
303 throw RipleyException(msg.str());
304 }
305 break;
306 default:
307 throw RipleyException(msg.str());
308 }
309 }
310
311 Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,
312 bool reducedColOrder) const
313 {
314 if (reducedRowOrder || reducedColOrder)
315 throw RipleyException("getPattern() not implemented for reduced order");
316
317 // connector
318 RankVector neighbour;
319 IndexVector offsetInShared(1,0);
320 IndexVector sendShared, recvShared;
321 const IndexVector faces=getNumFacesPerBoundary();
322 const index_t left = (m_offset0==0 ? 0 : 1);
323 const index_t bottom = (m_offset1==0 ? 0 : 1);
324 // corner node from bottom-left
325 if (faces[0] == 0 && faces[2] == 0) {
326 neighbour.push_back(m_mpiInfo->rank-m_NX-1);
327 offsetInShared.push_back(offsetInShared.back()+1);
328 sendShared.push_back(m_nodeId[m_N0+left]);
329 recvShared.push_back(m_nodeId[0]);
330 }
331 // bottom edge
332 if (faces[2] == 0) {
333 neighbour.push_back(m_mpiInfo->rank-m_NX);
334 offsetInShared.push_back(offsetInShared.back()+m_N0-left);
335 for (dim_t i=left; i<m_N0; i++) {
336 // easy case, we know the neighbour id's
337 sendShared.push_back(m_nodeId[m_N0+i]);
338 recvShared.push_back(m_nodeId[i]);
339 }
340 }
341 // corner node from bottom-right
342 if (faces[1] == 0 && faces[2] == 0) {
343 neighbour.push_back(m_mpiInfo->rank-m_NX+1);
344 const index_t N0=(neighbour.back()%m_NX == 0 ? m_N0 : m_N0-1);
345 const index_t N1=(neighbour.back()/m_NX == 0 ? m_N1 : m_N1-1);
346 const index_t first=m_nodeDistribution[neighbour.back()];
347 offsetInShared.push_back(offsetInShared.back()+1);
348 sendShared.push_back(m_nodeId[(bottom+1)*m_N0-1]);
349 recvShared.push_back(first+N0*(N1-1));
350 }
351 // left edge
352 if (faces[0] == 0) {
353 neighbour.push_back(m_mpiInfo->rank-1);
354 offsetInShared.push_back(offsetInShared.back()+m_N1-bottom);
355 for (dim_t i=bottom; i<m_N1; i++) {
356 // easy case, we know the neighbour id's
357 sendShared.push_back(m_nodeId[i*m_N0+1]);
358 recvShared.push_back(m_nodeId[i*m_N0]);
359 }
360 }
361 // right edge
362 if (faces[1] == 0) {
363 neighbour.push_back(m_mpiInfo->rank+1);
364 const index_t rightN0=(neighbour.back()%m_NX == 0 ? m_N0 : m_N0-1);
365 const index_t first=m_nodeDistribution[neighbour.back()];
366 offsetInShared.push_back(offsetInShared.back()+m_N1-bottom);
367 for (dim_t i=bottom; i<m_N1; i++) {
368 sendShared.push_back(m_nodeId[(i+1)*m_N0-1]);
369 recvShared.push_back(first+rightN0*(i-bottom));
370 }
371 }
372 // corner node from top-left
373 if (faces[0] == 0 && faces[3] == 0) {
374 neighbour.push_back(m_mpiInfo->rank+m_NX-1);
375 const index_t N0=(neighbour.back()%m_NX == 0 ? m_N0 : m_N0-1);
376 const index_t first=m_nodeDistribution[neighbour.back()];
377 offsetInShared.push_back(offsetInShared.back()+1);
378 sendShared.push_back(m_nodeId[m_N0*(m_N1-1)+left]);
379 recvShared.push_back(first+N0-1);
380 }
381 // top edge
382 if (faces[3] == 0) {
383 neighbour.push_back(m_mpiInfo->rank+m_NX);
384 const index_t first=m_nodeDistribution[neighbour.back()];
385 offsetInShared.push_back(offsetInShared.back()+m_N0-left);
386 for (dim_t i=left; i<m_N0; i++) {
387 sendShared.push_back(m_nodeId[m_N0*(m_N1-1)+i]);
388 recvShared.push_back(first+i-left);
389 }
390 }
391 // corner node from top-right
392 if (faces[1] == 0 && faces[3] == 0) {
393 neighbour.push_back(m_mpiInfo->rank+m_NX+1);
394 const index_t first=m_nodeDistribution[neighbour.back()];
395 offsetInShared.push_back(offsetInShared.back()+1);
396 sendShared.push_back(m_nodeId[m_N0*m_N1-1]);
397 recvShared.push_back(first);
398 }
399 const int numDOF=m_nodeDistribution[m_mpiInfo->rank+1]-m_nodeDistribution[m_mpiInfo->rank];
400 cout << "--- rcv_shcomp ---" << endl;
401 cout << "numDOF=" << numDOF << ", numNeighbors=" << neighbour.size() << endl;
402 for (size_t i=0; i<neighbour.size(); i++) {
403 cout << "neighbor[" << i << "]=" << neighbour[i]
404 << " offsetInShared[" << i+1 << "]=" << offsetInShared[i+1] << endl;
405 }
406 for (size_t i=0; i<recvShared.size(); i++) {
407 cout << "shared[" << i << "]=" << recvShared[i] << endl;
408 }
409 cout << "--- snd_shcomp ---" << endl;
410 for (size_t i=0; i<sendShared.size(); i++) {
411 cout << "shared[" << i << "]=" << sendShared[i] << endl;
412 }
413
414 Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
415 numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
416 &offsetInShared[0], 1, 0, m_mpiInfo);
417 Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(
418 numDOF, neighbour.size(), &neighbour[0], &recvShared[0],
419 &offsetInShared[0], 1, 0, m_mpiInfo);
420 Paso_Connector* connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);
421 Paso_SharedComponents_free(snd_shcomp);
422 Paso_SharedComponents_free(rcv_shcomp);
423
424 // create patterns
425 dim_t M, N;
426 IndexVector ptr(1,0);
427 IndexVector index;
428
429 // main pattern
430 for (index_t i=0; i<numDOF; i++) {
431 // always add the node itself
432 index.push_back(i);
433 int num=insertNeighbours(index, i);
434 ptr.push_back(ptr.back()+num+1);
435 }
436 M=N=ptr.size()-1;
437 // paso will manage the memory
438 index_t* indexC = MEMALLOC(index.size(),index_t);
439 index_t* ptrC = MEMALLOC(ptr.size(), index_t);
440 copy(index.begin(), index.end(), indexC);
441 copy(ptr.begin(), ptr.end(), ptrC);
442 Paso_Pattern *mainPattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,
443 M, N, ptrC, indexC);
444
445 cout << "--- main_pattern ---" << endl;
446 cout << "M=" << M << ", N=" << N << endl;
447 for (size_t i=0; i<ptr.size(); i++) {
448 cout << "ptr[" << i << "]=" << ptr[i] << endl;
449 }
450 for (size_t i=0; i<index.size(); i++) {
451 cout << "index[" << i << "]=" << index[i] << endl;
452 }
453
454 ptr.clear();
455 index.clear();
456
457 // column & row couple patterns
458 Paso_Pattern *colCouplePattern, *rowCouplePattern;
459 generateCouplePatterns(&colCouplePattern, &rowCouplePattern);
460
461 // allocate paso distribution
462 Paso_Distribution* distribution = Paso_Distribution_alloc(m_mpiInfo,
463 const_cast<index_t*>(&m_nodeDistribution[0]), 1, 0);
464
465 Paso_SystemMatrixPattern* pattern = Paso_SystemMatrixPattern_alloc(
466 MATRIX_FORMAT_DEFAULT, distribution, distribution,
467 mainPattern, colCouplePattern, rowCouplePattern,
468 connector, connector);
469 Paso_Pattern_free(mainPattern);
470 Paso_Pattern_free(colCouplePattern);
471 Paso_Pattern_free(rowCouplePattern);
472 Paso_Distribution_free(distribution);
473 return pattern;
474 }
475
476 void Rectangle::Print_Mesh_Info(const bool full) const
477 {
478 RipleyDomain::Print_Mesh_Info(full);
479 if (full) {
480 cout << " Id Coordinates" << endl;
481 cout.precision(15);
482 cout.setf(ios::scientific, ios::floatfield);
483 pair<double,double> xdx = getFirstCoordAndSpacing(0);
484 pair<double,double> ydy = getFirstCoordAndSpacing(1);
485 for (index_t i=0; i < getNumNodes(); i++) {
486 cout << " " << setw(5) << m_nodeId[i]
487 << " " << xdx.first+(i%m_N0)*xdx.second
488 << " " << ydy.first+(i/m_N0)*ydy.second << endl;
489 }
490 }
491 }
492
493 IndexVector Rectangle::getNumNodesPerDim() const
494 {
495 IndexVector ret;
496 ret.push_back(m_N0);
497 ret.push_back(m_N1);
498 return ret;
499 }
500
501 IndexVector Rectangle::getNumElementsPerDim() const
502 {
503 IndexVector ret;
504 ret.push_back(m_NE0);
505 ret.push_back(m_NE1);
506 return ret;
507 }
508
509 IndexVector Rectangle::getNumFacesPerBoundary() const
510 {
511 IndexVector ret(4, 0);
512 //left
513 if (m_offset0==0)
514 ret[0]=m_NE1;
515 //right
516 if (m_mpiInfo->rank%m_NX==m_NX-1)
517 ret[1]=m_NE1;
518 //bottom
519 if (m_offset1==0)
520 ret[2]=m_NE0;
521 //top
522 if (m_mpiInfo->rank/m_NX==m_NY-1)
523 ret[3]=m_NE0;
524 return ret;
525 }
526
527 pair<double,double> Rectangle::getFirstCoordAndSpacing(dim_t dim) const
528 {
529 if (dim==0) {
530 return pair<double,double>((m_l0*m_offset0)/m_gNE0, m_l0/m_gNE0);
531 } else if (dim==1) {
532 return pair<double,double>((m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);
533 }
534 throw RipleyException("getFirstCoordAndSpacing(): invalid argument");
535 }
536
537 //protected
538 dim_t Rectangle::getNumFaceElements() const
539 {
540 const IndexVector faces = getNumFacesPerBoundary();
541 dim_t n=0;
542 for (size_t i=0; i<faces.size(); i++)
543 n+=faces[i];
544 return n;
545 }
546
547 //protected
548 void Rectangle::assembleCoordinates(escript::Data& arg) const
549 {
550 escriptDataC x = arg.getDataC();
551 int numDim = m_numDim;
552 if (!isDataPointShapeEqual(&x, 1, &numDim))
553 throw RipleyException("setToX: Invalid Data object shape");
554 if (!numSamplesEqual(&x, 1, getNumNodes()))
555 throw RipleyException("setToX: Illegal number of samples in Data object");
556
557 pair<double,double> xdx = getFirstCoordAndSpacing(0);
558 pair<double,double> ydy = getFirstCoordAndSpacing(1);
559 arg.requireWrite();
560 #pragma omp parallel for
561 for (dim_t i1 = 0; i1 < m_N1; i1++) {
562 for (dim_t i0 = 0; i0 < m_N0; i0++) {
563 double* point = arg.getSampleDataRW(i0+m_N0*i1);
564 point[0] = xdx.first+i0*xdx.second;
565 point[1] = ydy.first+i1*ydy.second;
566 }
567 }
568 }
569
570 //private
571 void Rectangle::populateSampleIds()
572 {
573 // identifiers are ordered from left to right, bottom to top on each rank,
574 // except for the shared nodes which are owned by the rank below / to the
575 // left of the current rank
576
577 // build node distribution vector first.
578 // m_nodeDistribution[i] is the first node id on rank i, that is
579 // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes
580 m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
581 m_nodeDistribution[1]=getNumNodes();
582 for (dim_t k=1; k<m_mpiInfo->size-1; k++) {
583 const index_t x=k%m_NX;
584 const index_t y=k/m_NX;
585 index_t numNodes=getNumNodes();
586 if (x>0)
587 numNodes-=m_N1;
588 if (y>0)
589 numNodes-=m_N0;
590 if (x>0 && y>0)
591 numNodes++; // subtracted corner twice -> fix that
592 m_nodeDistribution[k+1]=m_nodeDistribution[k]+numNodes;
593 }
594 m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
595
596 m_nodeId.resize(getNumNodes());
597
598 // the bottom row and left column are not owned by this rank so the
599 // identifiers need to be computed accordingly
600 const index_t left = (m_offset0==0 ? 0 : 1);
601 const index_t bottom = (m_offset1==0 ? 0 : 1);
602 if (left>0) {
603 const int neighbour=m_mpiInfo->rank-1;
604 const index_t leftN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
605 #pragma omp parallel for
606 for (dim_t i1=bottom; i1<m_N1; i1++) {
607 m_nodeId[i1*m_N0]=m_nodeDistribution[neighbour]
608 + (i1-bottom+1)*leftN0-1;
609 }
610 }
611 if (bottom>0) {
612 const int neighbour=m_mpiInfo->rank-m_NX;
613 const index_t bottomN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
614 const index_t bottomN1=(neighbour/m_NX == 0 ? m_N1 : m_N1-1);
615 #pragma omp parallel for
616 for (dim_t i0=left; i0<m_N0; i0++) {
617 m_nodeId[i0]=m_nodeDistribution[neighbour]
618 + (bottomN1-1)*bottomN0 + i0 - left;
619 }
620 }
621 if (left>0 && bottom>0) {
622 const int neighbour=m_mpiInfo->rank-m_NX-1;
623 const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
624 const index_t N1=(neighbour/m_NX == 0 ? m_N1 : m_N1-1);
625 m_nodeId[0]=m_nodeDistribution[neighbour]+N0*N1-1;
626 }
627
628 // the rest of the id's are contiguous
629 const index_t firstId=m_nodeDistribution[m_mpiInfo->rank];
630 #pragma omp parallel for
631 for (dim_t i1=bottom; i1<m_N1; i1++) {
632 for (dim_t i0=left; i0<m_N0; i0++) {
633 m_nodeId[i0+i1*m_N0] = firstId+i0-left+(i1-bottom)*(m_N0-left);
634 }
635 }
636
637 // elements
638 m_elementId.resize(getNumElements());
639 #pragma omp parallel for
640 for (dim_t k=0; k<getNumElements(); k++) {
641 m_elementId[k]=k;
642 }
643
644 // face elements
645 m_faceId.resize(getNumFaceElements());
646 #pragma omp parallel for
647 for (dim_t k=0; k<getNumFaceElements(); k++) {
648 m_faceId[k]=k;
649 }
650 }
651
652 //private
653 int Rectangle::insertNeighbours(IndexVector& index, index_t node) const
654 {
655 const dim_t myN0 = (m_offset0==0 ? m_N0 : m_N0-1);
656 const dim_t myN1 = (m_offset1==0 ? m_N1 : m_N1-1);
657 const int x=node%myN0;
658 const int y=node/myN0;
659 int num=0;
660 if (y>0) {
661 if (x>0) {
662 // bottom-left
663 index.push_back(node-myN0-1);
664 num++;
665 }
666 // bottom
667 index.push_back(node-myN0);
668 num++;
669 if (x<myN0-1) {
670 // bottom-right
671 index.push_back(node-myN0+1);
672 num++;
673 }
674 }
675 if (x<myN0-1) {
676 // right
677 index.push_back(node+1);
678 num++;
679 if (y<myN1-1) {
680 // top-right
681 index.push_back(node+myN0+1);
682 num++;
683 }
684 }
685 if (y<myN1-1) {
686 // top
687 index.push_back(node+myN0);
688 num++;
689 if (x>0) {
690 // top-left
691 index.push_back(node+myN0-1);
692 num++;
693 }
694 }
695 if (x>0) {
696 // left
697 index.push_back(node-1);
698 num++;
699 }
700
701 return num;
702 }
703
704 //private
705 void Rectangle::generateCouplePatterns(Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const
706 {
707 IndexVector ptr(1,0);
708 IndexVector index;
709 const dim_t myN0 = (m_offset0==0 ? m_N0 : m_N0-1);
710 const dim_t myN1 = (m_offset1==0 ? m_N1 : m_N1-1);
711 const IndexVector faces=getNumFacesPerBoundary();
712
713 // bottom edge
714 dim_t n=0;
715 if (faces[0] == 0) {
716 index.push_back(2*(myN0+myN1+1));
717 index.push_back(2*(myN0+myN1+1)+1);
718 n+=2;
719 if (faces[2] == 0) {
720 index.push_back(0);
721 index.push_back(1);
722 index.push_back(2);
723 n+=3;
724 }
725 } else if (faces[2] == 0) {
726 index.push_back(1);
727 index.push_back(2);
728 n+=2;
729 }
730 // n=neighbours of bottom-left corner node
731 ptr.push_back(ptr.back()+n);
732 n=0;
733 if (faces[2] == 0) {
734 for (dim_t i=1; i<myN0-1; i++) {
735 index.push_back(i);
736 index.push_back(i+1);
737 index.push_back(i+2);
738 ptr.push_back(ptr.back()+3);
739 }
740 index.push_back(myN0-1);
741 index.push_back(myN0);
742 n+=2;
743 if (faces[1] == 0) {
744 index.push_back(myN0+1);
745 index.push_back(myN0+2);
746 index.push_back(myN0+3);
747 n+=3;
748 }
749 } else {
750 for (dim_t i=1; i<myN0-1; i++) {
751 ptr.push_back(ptr.back());
752 }
753 if (faces[1] == 0) {
754 index.push_back(myN0+2);
755 index.push_back(myN0+3);
756 n+=2;
757 }
758 }
759 // n=neighbours of bottom-right corner node
760 ptr.push_back(ptr.back()+n);
761
762 // 2nd row to 2nd last row
763 for (dim_t i=1; i<myN1-1; i++) {
764 // left edge
765 if (faces[0] == 0) {
766 index.push_back(2*(myN0+myN1+2)-i);
767 index.push_back(2*(myN0+myN1+2)-i-1);
768 index.push_back(2*(myN0+myN1+2)-i-2);
769 ptr.push_back(ptr.back()+3);
770 } else {
771 ptr.push_back(ptr.back());
772 }
773 for (dim_t j=1; j<myN0-1; j++) {
774 ptr.push_back(ptr.back());
775 }
776 // right edge
777 if (faces[1] == 0) {
778 index.push_back(myN0+i+1);
779 index.push_back(myN0+i+2);
780 index.push_back(myN0+i+3);
781 ptr.push_back(ptr.back()+3);
782 } else {
783 ptr.push_back(ptr.back());
784 }
785 }
786
787 // top edge
788 n=0;
789 if (faces[0] == 0) {
790 index.push_back(2*myN0+myN1+5);
791 index.push_back(2*myN0+myN1+4);
792 n+=2;
793 if (faces[3] == 0) {
794 index.push_back(2*myN0+myN1+3);
795 index.push_back(2*myN0+myN1+2);
796 index.push_back(2*myN0+myN1+1);
797 n+=3;
798 }
799 } else if (faces[3] == 0) {
800 index.push_back(2*myN0+myN1+2);
801 index.push_back(2*myN0+myN1+1);
802 n+=2;
803 }
804 // n=neighbours of top-left corner node
805 ptr.push_back(ptr.back()+n);
806 n=0;
807 if (faces[3] == 0) {
808 for (dim_t i=1; i<myN0-1; i++) {
809 index.push_back(2*myN0+myN1+i+1);
810 index.push_back(2*myN0+myN1+i);
811 index.push_back(2*myN0+myN1+i-1);
812 ptr.push_back(ptr.back()+3);
813 }
814 index.push_back(myN0+myN1+4);
815 index.push_back(myN0+myN1+3);
816 n+=2;
817 if (faces[1] == 0) {
818 index.push_back(myN0+myN1+2);
819 index.push_back(myN0+myN1+1);
820 index.push_back(myN0+myN1);
821 n+=3;
822 }
823 } else {
824 for (dim_t i=1; i<myN0-1; i++) {
825 ptr.push_back(ptr.back());
826 }
827 if (faces[1] == 0) {
828 index.push_back(myN0+myN1+1);
829 index.push_back(myN0+myN1);
830 n+=2;
831 }
832 }
833 // n=neighbours of top-right corner node
834 ptr.push_back(ptr.back()+n);
835
836 dim_t M=ptr.size()-1;
837 map<index_t,index_t> idMap;
838 dim_t N=0;
839 for (IndexVector::iterator it=index.begin(); it!=index.end(); it++) {
840 if (idMap.find(*it)==idMap.end()) {
841 idMap[*it]=N;
842 *it=N++;
843 } else {
844 *it=idMap[*it];
845 }
846 }
847
848 cout << "--- colCouple_pattern ---" << endl;
849 cout << "M=" << M << ", N=" << N << endl;
850 for (size_t i=0; i<ptr.size(); i++) {
851 cout << "ptr[" << i << "]=" << ptr[i] << endl;
852 }
853 for (size_t i=0; i<index.size(); i++) {
854 cout << "index[" << i << "]=" << index[i] << endl;
855 }
856
857 // now build the row couple pattern
858 IndexVector ptr2(1,0);
859 IndexVector index2;
860 for (dim_t id=0; id<N; id++) {
861 n=0;
862 for (dim_t i=0; i<M; i++) {
863 for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {
864 if (index[j]==id) {
865 index2.push_back(i);
866 n++;
867 break;
868 }
869 }
870 }
871 ptr2.push_back(ptr2.back()+n);
872 }
873
874 cout << "--- rowCouple_pattern ---" << endl;
875 cout << "M=" << N << ", N=" << M << endl;
876 for (size_t i=0; i<ptr2.size(); i++) {
877 cout << "ptr[" << i << "]=" << ptr2[i] << endl;
878 }
879 for (size_t i=0; i<index2.size(); i++) {
880 cout << "index[" << i << "]=" << index2[i] << endl;
881 }
882
883 // paso will manage the memory
884 index_t* indexC = MEMALLOC(index.size(), index_t);
885 index_t* ptrC = MEMALLOC(ptr.size(), index_t);
886 copy(index.begin(), index.end(), indexC);
887 copy(ptr.begin(), ptr.end(), ptrC);
888 *colPattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);
889
890 // paso will manage the memory
891 indexC = MEMALLOC(index2.size(), index_t);
892 ptrC = MEMALLOC(ptr2.size(), index_t);
893 copy(index2.begin(), index2.end(), indexC);
894 copy(ptr2.begin(), ptr2.end(), ptrC);
895 *rowPattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, N, M, ptrC, indexC);
896 }
897
898 } // end of namespace ripley
899

  ViewVC Help
Powered by ViewVC 1.1.26