/[escript]/branches/ripleygmg_from_3668/ripley/src/Rectangle.cpp
ViewVC logotype

Annotation of /branches/ripleygmg_from_3668/ripley/src/Rectangle.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3699 - (hide annotations)
Thu Dec 1 22:59:42 2011 UTC (8 years ago) by caltinay
File size: 29713 byte(s)
getPattern() for Rectangle. Couple node numbering is likely to change once
we get to the PDE assembly.

1 caltinay 3691
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 caltinay 3697 pair<double,double> xdx = getFirstCoordAndSpacing(0);
142     pair<double,double> ydy = getFirstCoordAndSpacing(1);
143 caltinay 3691 #pragma omp parallel
144     {
145     #pragma omp for
146     for (dim_t i0 = 0; i0 < m_N0; i0++) {
147 caltinay 3697 coords[0][i0]=xdx.first+i0*xdx.second;
148 caltinay 3691 }
149     #pragma omp for
150     for (dim_t i1 = 0; i1 < m_N1; i1++) {
151 caltinay 3697 coords[1][i1]=ydy.first+i1*ydy.second;
152 caltinay 3691 }
153     }
154 caltinay 3697 IndexVector dims = getNumNodesPerDim();
155    
156     // write mesh
157     DBPutQuadmesh(dbfile, "mesh", NULL, coords, &dims[0], 2, DB_DOUBLE,
158 caltinay 3691 DB_COLLINEAR, NULL);
159    
160 caltinay 3697 // 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 caltinay 3691 if (m_mpiInfo->rank == 0) {
171     vector<string> tempstrings;
172 caltinay 3697 vector<char*> names;
173 caltinay 3691 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 caltinay 3697 names.push_back((char*)tempstrings.back().c_str());
178 caltinay 3691 }
179 caltinay 3697 vector<int> types(m_mpiInfo->size, DB_QUAD_RECT);
180 caltinay 3691 DBSetDir(dbfile, "/");
181 caltinay 3697 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 caltinay 3691 }
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 caltinay 3697 const int* Rectangle::borrowSampleReferenceIDs(int fsType) const
221 caltinay 3691 {
222 caltinay 3697 switch (fsType) {
223 caltinay 3691 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 caltinay 3697 << fsType;
236 caltinay 3691 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 caltinay 3699 const index_t myFirst=m_nodeDistribution[m_mpiInfo->rank];
244     const index_t myLast=m_nodeDistribution[m_mpiInfo->rank+1]-1;
245 caltinay 3691 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 caltinay 3697 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 caltinay 3699 case DegreesOfFreedom:
273     target=in;
274     break;
275    
276 caltinay 3697 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 caltinay 3699 break;
306 caltinay 3697 default:
307     throw RipleyException(msg.str());
308     }
309 caltinay 3691 }
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 caltinay 3699 // 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 caltinay 3691
414     Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
415 caltinay 3699 numDOF, neighbour.size(), &neighbour[0], &sendShared[0],
416 caltinay 3697 &offsetInShared[0], 1, 0, m_mpiInfo);
417 caltinay 3691 Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(
418 caltinay 3699 numDOF, neighbour.size(), &neighbour[0], &recvShared[0],
419 caltinay 3697 &offsetInShared[0], 1, 0, m_mpiInfo);
420 caltinay 3691 Paso_Connector* connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);
421 caltinay 3699 Paso_SharedComponents_free(snd_shcomp);
422     Paso_SharedComponents_free(rcv_shcomp);
423 caltinay 3691
424     // create patterns
425 caltinay 3699 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 caltinay 3691 Paso_Pattern *mainPattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,
443 caltinay 3699 M, N, ptrC, indexC);
444 caltinay 3691
445 caltinay 3699 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 caltinay 3691 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 caltinay 3697 return pattern;
474 caltinay 3691 }
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 caltinay 3697 pair<double,double> xdx = getFirstCoordAndSpacing(0);
484     pair<double,double> ydy = getFirstCoordAndSpacing(1);
485 caltinay 3691 for (index_t i=0; i < getNumNodes(); i++) {
486     cout << " " << setw(5) << m_nodeId[i]
487 caltinay 3697 << " " << xdx.first+(i%m_N0)*xdx.second
488     << " " << ydy.first+(i/m_N0)*ydy.second << endl;
489 caltinay 3691 }
490     }
491     }
492    
493 caltinay 3697 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 caltinay 3691 //protected
538     dim_t Rectangle::getNumFaceElements() const
539     {
540 caltinay 3699 const IndexVector faces = getNumFacesPerBoundary();
541 caltinay 3691 dim_t n=0;
542 caltinay 3699 for (size_t i=0; i<faces.size(); i++)
543     n+=faces[i];
544 caltinay 3691 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 caltinay 3697 pair<double,double> xdx = getFirstCoordAndSpacing(0);
558     pair<double,double> ydy = getFirstCoordAndSpacing(1);
559 caltinay 3691 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 caltinay 3697 point[0] = xdx.first+i0*xdx.second;
565     point[1] = ydy.first+i1*ydy.second;
566 caltinay 3691 }
567     }
568     }
569    
570     //private
571     void Rectangle::populateSampleIds()
572     {
573 caltinay 3697 // 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 caltinay 3691 m_nodeId.resize(getNumNodes());
597 caltinay 3697
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 caltinay 3691 #pragma omp parallel for
606 caltinay 3697 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 caltinay 3699 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 caltinay 3697 }
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 caltinay 3691 }
651    
652 caltinay 3699 //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 caltinay 3691 } // end of namespace ripley
899    

  ViewVC Help
Powered by ViewVC 1.1.26