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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3698 - (hide annotations)
Tue Nov 29 00:47:29 2011 UTC (7 years, 4 months ago) by caltinay
File size: 18928 byte(s)
[x] Rectangle...
[x] Brick node id's and weipa compatibility


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/Brick.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     Brick::Brick(int n0, int n1, int n2, double l0, double l1, double l2, int d0,
33     int d1, int d2) :
34     RipleyDomain(3),
35     m_gNE0(n0),
36     m_gNE1(n1),
37     m_gNE2(n2),
38     m_l0(l0),
39     m_l1(l1),
40     m_l2(l2),
41     m_NX(d0),
42     m_NY(d1),
43     m_NZ(d2)
44     {
45     // ensure number of subdivisions is valid and nodes can be distributed
46     // among number of ranks
47     if (m_NX*m_NY*m_NZ != m_mpiInfo->size)
48     throw RipleyException("Invalid number of spatial subdivisions");
49    
50     if (n0%m_NX > 0 || n1%m_NY > 0 || n2%m_NZ > 0)
51     throw RipleyException("Number of elements must be separable into number of ranks in each dimension");
52    
53     // local number of elements
54     m_NE0 = n0/m_NX;
55     m_NE1 = n1/m_NY;
56     m_NE2 = n2/m_NZ;
57     // local number of nodes (not necessarily owned)
58     m_N0 = m_NE0+1;
59     m_N1 = m_NE1+1;
60     m_N2 = m_NE2+1;
61     // bottom-left-front node is at (offset0,offset1,offset2) in global mesh
62     m_offset0 = m_NE0*(m_mpiInfo->rank%m_NX);
63     m_offset1 = m_NE1*(m_mpiInfo->rank%(m_NX*m_NY)/m_NX);
64     m_offset2 = m_NE2*(m_mpiInfo->rank/(m_NX*m_NY));
65     populateSampleIds();
66     }
67    
68    
69     Brick::~Brick()
70     {
71     }
72    
73     string Brick::getDescription() const
74     {
75     return "ripley::Brick";
76     }
77    
78     bool Brick::operator==(const AbstractDomain& other) const
79     {
80     if (dynamic_cast<const Brick*>(&other))
81     return this==&other;
82    
83     return false;
84     }
85    
86     void Brick::dump(const string& fileName) const
87     {
88     #if USE_SILO
89     string fn(fileName);
90     if (fileName.length() < 6 || fileName.compare(fileName.length()-5, 5, ".silo") != 0) {
91     fn+=".silo";
92     }
93    
94     const int NUM_SILO_FILES = 1;
95     const char* blockDirFmt = "/block%04d";
96     int driver=DB_HDF5;
97     string siloPath;
98     DBfile* dbfile = NULL;
99    
100     #ifdef ESYS_MPI
101     PMPIO_baton_t* baton = NULL;
102     #endif
103    
104     if (m_mpiInfo->size > 1) {
105     #ifdef ESYS_MPI
106     baton = PMPIO_Init(NUM_SILO_FILES, PMPIO_WRITE, m_mpiInfo->comm,
107     0x1337, PMPIO_DefaultCreate, PMPIO_DefaultOpen,
108     PMPIO_DefaultClose, (void*)&driver);
109     // try the fallback driver in case of error
110     if (!baton && driver != DB_PDB) {
111     driver = DB_PDB;
112     baton = PMPIO_Init(NUM_SILO_FILES, PMPIO_WRITE, m_mpiInfo->comm,
113     0x1338, PMPIO_DefaultCreate, PMPIO_DefaultOpen,
114     PMPIO_DefaultClose, (void*)&driver);
115     }
116     if (baton) {
117     char str[64];
118     snprintf(str, 64, blockDirFmt, PMPIO_RankInGroup(baton, m_mpiInfo->rank));
119     siloPath = str;
120     dbfile = (DBfile*) PMPIO_WaitForBaton(baton, fn.c_str(), siloPath.c_str());
121     }
122     #endif
123     } else {
124     dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,
125     getDescription().c_str(), driver);
126     // try the fallback driver in case of error
127     if (!dbfile && driver != DB_PDB) {
128     driver = DB_PDB;
129     dbfile = DBCreate(fn.c_str(), DB_CLOBBER, DB_LOCAL,
130     getDescription().c_str(), driver);
131     }
132     }
133    
134     if (!dbfile)
135     throw RipleyException("dump: Could not create Silo file");
136    
137     /*
138     if (driver==DB_HDF5) {
139     // gzip level 1 already provides good compression with minimal
140     // performance penalty. Some tests showed that gzip levels >3 performed
141     // rather badly on escript data both in terms of time and space
142     DBSetCompression("ERRMODE=FALLBACK METHOD=GZIP LEVEL=1");
143     }
144     */
145    
146     boost::scoped_ptr<double> x(new double[m_N0]);
147     boost::scoped_ptr<double> y(new double[m_N1]);
148     boost::scoped_ptr<double> z(new double[m_N2]);
149     double* coords[3] = { x.get(), y.get(), z.get() };
150 caltinay 3698 pair<double,double> xdx = getFirstCoordAndSpacing(0);
151     pair<double,double> ydy = getFirstCoordAndSpacing(1);
152     pair<double,double> zdz = getFirstCoordAndSpacing(2);
153 caltinay 3691 #pragma omp parallel
154     {
155     #pragma omp for
156     for (dim_t i0 = 0; i0 < m_N0; i0++) {
157 caltinay 3698 coords[0][i0]=xdx.first+i0*xdx.second;
158 caltinay 3691 }
159     #pragma omp for
160     for (dim_t i1 = 0; i1 < m_N1; i1++) {
161 caltinay 3698 coords[1][i1]=ydy.first+i1*ydy.second;
162 caltinay 3691 }
163     #pragma omp for
164     for (dim_t i2 = 0; i2 < m_N2; i2++) {
165 caltinay 3698 coords[2][i2]=zdz.first+i2*zdz.second;
166 caltinay 3691 }
167     }
168 caltinay 3698 IndexVector dims = getNumNodesPerDim();
169     DBPutQuadmesh(dbfile, "mesh", NULL, coords, &dims[0], 3, DB_DOUBLE,
170 caltinay 3691 DB_COLLINEAR, NULL);
171    
172 caltinay 3698 DBPutQuadvar1(dbfile, "nodeId", "mesh", (void*)&m_nodeId[0], &dims[0], 3,
173     NULL, 0, DB_INT, DB_NODECENT, NULL);
174 caltinay 3691
175 caltinay 3698 // write element ids
176     dims = getNumElementsPerDim();
177     DBPutQuadvar1(dbfile, "elementId", "mesh", (void*)&m_elementId[0],
178     &dims[0], 3, NULL, 0, DB_INT, DB_ZONECENT, NULL);
179    
180     // rank 0 writes multimesh and multivar
181 caltinay 3691 if (m_mpiInfo->rank == 0) {
182     vector<string> tempstrings;
183     vector<char*> names;
184     for (dim_t i=0; i<m_mpiInfo->size; i++) {
185     stringstream path;
186     path << "/block" << setw(4) << setfill('0') << right << i << "/mesh";
187     tempstrings.push_back(path.str());
188     names.push_back((char*)tempstrings.back().c_str());
189     }
190     vector<int> types(m_mpiInfo->size, DB_QUAD_RECT);
191     DBSetDir(dbfile, "/");
192     DBPutMultimesh(dbfile, "multimesh", 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 << "/nodeId";
199     tempstrings.push_back(path.str());
200     names.push_back((char*)tempstrings.back().c_str());
201     }
202     types.assign(m_mpiInfo->size, DB_QUADVAR);
203     DBPutMultivar(dbfile, "nodeId", m_mpiInfo->size, &names[0],
204     &types[0], NULL);
205 caltinay 3698 tempstrings.clear();
206     names.clear();
207     for (dim_t i=0; i<m_mpiInfo->size; i++) {
208     stringstream path;
209     path << "/block" << setw(4) << setfill('0') << right << i << "/elementId";
210     tempstrings.push_back(path.str());
211     names.push_back((char*)tempstrings.back().c_str());
212     }
213     DBPutMultivar(dbfile, "elementId", m_mpiInfo->size, &names[0],
214     &types[0], NULL);
215 caltinay 3691 }
216    
217     if (m_mpiInfo->size > 1) {
218     #ifdef ESYS_MPI
219     PMPIO_HandOffBaton(baton, dbfile);
220     PMPIO_Finish(baton);
221     #endif
222     } else {
223     DBClose(dbfile);
224     }
225    
226     #else // USE_SILO
227     throw RipleyException("dump(): no Silo support");
228     #endif
229     }
230    
231     const int* Brick::borrowSampleReferenceIDs(int fsType) const
232     {
233 caltinay 3697 switch (fsType) {
234     case Nodes:
235     return &m_nodeId[0];
236     case Elements:
237     return &m_elementId[0];
238     case FaceElements:
239     return &m_faceId[0];
240     default:
241     break;
242     }
243 caltinay 3691
244 caltinay 3697 stringstream msg;
245     msg << "borrowSampleReferenceIDs() not implemented for function space type "
246     << fsType;
247     throw RipleyException(msg.str());
248 caltinay 3691 }
249    
250     bool Brick::ownSample(int fsCode, index_t id) const
251     {
252     #ifdef ESYS_MPI
253     if (fsCode == Nodes) {
254 caltinay 3698 const index_t myFirst=m_nodeDistribution[m_mpiInfo->rank];
255     const index_t myLast=m_nodeDistribution[m_mpiInfo->rank+1]-1;
256 caltinay 3691 return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);
257     } else
258     throw RipleyException("ownSample() only implemented for Nodes");
259     #else
260     return true;
261     #endif
262     }
263    
264     void Brick::interpolateOnDomain(escript::Data& target,
265     const escript::Data& in) const
266     {
267     const Brick& inDomain=dynamic_cast<const Brick&>(*(in.getFunctionSpace().getDomain()));
268     const Brick& targetDomain=dynamic_cast<const Brick&>(*(target.getFunctionSpace().getDomain()));
269     if (inDomain != *this)
270     throw RipleyException("Illegal domain of interpolant");
271     if (targetDomain != *this)
272     throw RipleyException("Illegal domain of interpolation target");
273    
274     throw RipleyException("interpolateOnDomain() not implemented");
275     }
276    
277     Paso_SystemMatrixPattern* Brick::getPattern(bool reducedRowOrder,
278     bool reducedColOrder) const
279     {
280     if (reducedRowOrder || reducedColOrder)
281     throw RipleyException("getPattern() not implemented for reduced order");
282    
283     throw RipleyException("getPattern() not implemented");
284     }
285    
286     void Brick::Print_Mesh_Info(const bool full) const
287     {
288     RipleyDomain::Print_Mesh_Info(full);
289     if (full) {
290     cout << " Id Coordinates" << endl;
291     cout.precision(15);
292     cout.setf(ios::scientific, ios::floatfield);
293 caltinay 3698 pair<double,double> xdx = getFirstCoordAndSpacing(0);
294     pair<double,double> ydy = getFirstCoordAndSpacing(1);
295     pair<double,double> zdz = getFirstCoordAndSpacing(2);
296 caltinay 3691 for (index_t i=0; i < getNumNodes(); i++) {
297     cout << " " << setw(5) << m_nodeId[i]
298 caltinay 3698 << " " << xdx.first+(i%m_N0)*xdx.second
299     << " " << ydy.first+(i%(m_N0*m_N1)/m_N0)*ydy.second
300     << " " << zdz.first+(i/(m_N0*m_N1))*zdz.second << endl;
301 caltinay 3691 }
302     }
303     }
304    
305 caltinay 3698 IndexVector Brick::getNumNodesPerDim() const
306     {
307     IndexVector ret;
308     ret.push_back(m_N0);
309     ret.push_back(m_N1);
310     ret.push_back(m_N2);
311     return ret;
312     }
313    
314     IndexVector Brick::getNumElementsPerDim() const
315     {
316     IndexVector ret;
317     ret.push_back(m_NE0);
318     ret.push_back(m_NE1);
319     ret.push_back(m_NE2);
320     return ret;
321     }
322    
323     IndexVector Brick::getNumFacesPerBoundary() const
324     {
325     IndexVector ret(6, 0);
326     //left
327     if (m_offset0==0)
328     ret[0]=m_NE1*m_NE2;
329     //right
330     if (m_mpiInfo->rank%m_NX==m_NX-1)
331     ret[1]=m_NE1*m_NE2;
332     //bottom
333     if (m_offset1==0)
334     ret[2]=m_NE0*m_NE2;
335     //top
336     if (m_mpiInfo->rank%(m_NX*m_NY)/m_NX==m_NY-1)
337     ret[3]=m_NE0*m_NE2;
338     //front
339     if (m_offset2==0)
340     ret[4]=m_NE0*m_NE1;
341     //back
342     if (m_mpiInfo->rank/(m_NX*m_NY)==m_NZ-1)
343     ret[5]=m_NE0*m_NE1;
344     return ret;
345     }
346    
347     pair<double,double> Brick::getFirstCoordAndSpacing(dim_t dim) const
348     {
349     if (dim==0)
350     return pair<double,double>((m_l0*m_offset0)/m_gNE0, m_l0/m_gNE0);
351     else if (dim==1)
352     return pair<double,double>((m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);
353     else if (dim==2)
354     return pair<double,double>((m_l2*m_offset2)/m_gNE2, m_l2/m_gNE2);
355    
356     throw RipleyException("getFirstCoordAndSpacing(): invalid argument");
357     }
358    
359    
360 caltinay 3691 //protected
361     dim_t Brick::getNumFaceElements() const
362     {
363     dim_t n=0;
364     //left
365     if (m_offset0==0)
366     n+=m_NE1*m_NE2;
367     //right
368     if (m_mpiInfo->rank%m_NX==m_NX-1)
369     n+=m_NE1*m_NE2;
370     //bottom
371     if (m_offset1==0)
372     n+=m_NE0*m_NE2;
373     //top
374     if (m_mpiInfo->rank%(m_NX*m_NY)/m_NX==m_NY-1)
375     n+=m_NE0*m_NE2;
376     //front
377     if (m_offset2==0)
378     n+=m_NE0*m_NE1;
379     //back
380     if (m_mpiInfo->rank/(m_NX*m_NY)==m_NZ-1)
381     n+=m_NE0*m_NE1;
382    
383     return n;
384     }
385    
386     //protected
387     void Brick::assembleCoordinates(escript::Data& arg) const
388     {
389     escriptDataC x = arg.getDataC();
390     int numDim = m_numDim;
391     if (!isDataPointShapeEqual(&x, 1, &numDim))
392     throw RipleyException("setToX: Invalid Data object shape");
393     if (!numSamplesEqual(&x, 1, getNumNodes()))
394     throw RipleyException("setToX: Illegal number of samples in Data object");
395    
396 caltinay 3698 pair<double,double> xdx = getFirstCoordAndSpacing(0);
397     pair<double,double> ydy = getFirstCoordAndSpacing(1);
398     pair<double,double> zdz = getFirstCoordAndSpacing(2);
399 caltinay 3691 arg.requireWrite();
400     #pragma omp parallel for
401     for (dim_t i2 = 0; i2 < m_N2; i2++) {
402     for (dim_t i1 = 0; i1 < m_N1; i1++) {
403     for (dim_t i0 = 0; i0 < m_N0; i0++) {
404     double* point = arg.getSampleDataRW(i0+m_N0*i1+m_N0*m_N1*i2);
405 caltinay 3698 point[0] = xdx.first+i0*xdx.second;
406     point[1] = ydy.first+i1*ydy.second;
407     point[2] = zdz.first+i2*zdz.second;
408 caltinay 3691 }
409     }
410     }
411     }
412    
413     //private
414     void Brick::populateSampleIds()
415     {
416 caltinay 3697 // identifiers are ordered from left to right, bottom to top, front to back
417     // on each rank, except for the shared nodes which are owned by the rank
418     // below / to the left / to the front of the current rank
419 caltinay 3698
420     // build node distribution vector first.
421     // m_nodeDistribution[i] is the first node id on rank i, that is
422     // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes
423     m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
424     m_nodeDistribution[1]=getNumNodes();
425     for (dim_t k=1; k<m_mpiInfo->size-1; k++) {
426     const index_t x = k%m_NX;
427     const index_t y = k%(m_NX*m_NY)/m_NX;
428     const index_t z = k/(m_NX*m_NY);
429     index_t numNodes=getNumNodes();
430     if (x>0)
431     numNodes-=m_N1*m_N2;
432     if (y>0)
433     numNodes-=m_N0*m_N2;
434     if (z>0)
435     numNodes-=m_N0*m_N1;
436     // if an edge was subtracted twice add it back
437     if (x>0 && y>0)
438     numNodes+=m_N2;
439     if (x>0 && z>0)
440     numNodes+=m_N1;
441     if (y>0 && z>0)
442     numNodes+=m_N0;
443     // the corner node was removed 3x and added back 3x, so subtract it
444     if (x>0 && y>0 && z>0)
445     numNodes--;
446     m_nodeDistribution[k+1]=m_nodeDistribution[k]+numNodes;
447     }
448     m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
449    
450 caltinay 3691 m_nodeId.resize(getNumNodes());
451 caltinay 3698
452     // the bottom, left and front planes are not owned by this rank so the
453     // identifiers need to be computed accordingly
454     const index_t left = (m_offset0==0 ? 0 : 1);
455     const index_t bottom = (m_offset1==0 ? 0 : 1);
456     const index_t front = (m_offset2==0 ? 0 : 1);
457    
458     // case 1: all nodes on left plane are owned by rank on the left
459     if (left>0) {
460     const int neighbour=m_mpiInfo->rank-1;
461     const index_t leftN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
462     const index_t leftN1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
463 caltinay 3691 #pragma omp parallel for
464 caltinay 3698 for (dim_t i2=front; i2<m_N2; i2++) {
465     for (dim_t i1=bottom; i1<m_N1; i1++) {
466     m_nodeId[i1*m_N0+i2*m_N0*m_N1]=m_nodeDistribution[neighbour]
467     + (i1-bottom+1)*leftN0
468     + (i2-front)*leftN0*leftN1 - 1;
469     }
470     }
471 caltinay 3691 }
472 caltinay 3698 // case 2: all nodes on bottom plane are owned by rank below
473     if (bottom>0) {
474     const int neighbour=m_mpiInfo->rank-m_NX;
475     const index_t bottomN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
476     const index_t bottomN1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
477     #pragma omp parallel for
478     for (dim_t i2=front; i2<m_N2; i2++) {
479     for (dim_t i0=left; i0<m_N0; i0++) {
480     m_nodeId[i0+i2*m_N0*m_N1]=m_nodeDistribution[neighbour]
481     + bottomN0*(bottomN1-1)
482     + (i2-front)*bottomN0*bottomN1 + i0-left;
483     }
484     }
485     }
486     // case 3: all nodes on front plane are owned by rank in front
487     if (front>0) {
488     const int neighbour=m_mpiInfo->rank-m_NX*m_NY;
489     const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
490     const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
491     const index_t N2=(neighbour/(m_NX*m_NY)==0 ? m_N2 : m_N2-1);
492     #pragma omp parallel for
493     for (dim_t i1=bottom; i1<m_N1; i1++) {
494     for (dim_t i0=left; i0<m_N0; i0++) {
495     m_nodeId[i0+i1*m_N0]=m_nodeDistribution[neighbour]
496     + N0*N1*(N2-1)+(i1-bottom)*N0 + i0-left;
497     }
498     }
499     }
500     // case 4: nodes on front bottom edge are owned by the corresponding rank
501     if (front>0 && bottom>0) {
502     const int neighbour=m_mpiInfo->rank-m_NX*(m_NY+1);
503     const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
504     const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
505     const index_t N2=(neighbour/(m_NX*m_NY)==0 ? m_N2 : m_N2-1);
506     #pragma omp parallel for
507     for (dim_t i0=left; i0<m_N0; i0++) {
508     m_nodeId[i0]=m_nodeDistribution[neighbour]
509     + N0*N1*(N2-1)+(N1-1)*N0 + i0-left;
510     }
511     }
512     // case 5: nodes on left bottom edge are owned by the corresponding rank
513     if (left>0 && bottom>0) {
514     const int neighbour=m_mpiInfo->rank-m_NX-1;
515     const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
516     const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
517     #pragma omp parallel for
518     for (dim_t i2=front; i2<m_N2; i2++) {
519     m_nodeId[i2*m_N0*m_N1]=m_nodeDistribution[neighbour]
520     + (1+i2-front)*N0*N1-1;
521     }
522     }
523     // case 6: nodes on left front edge are owned by the corresponding rank
524     if (left>0 && front>0) {
525     const int neighbour=m_mpiInfo->rank-m_NX*m_NY-1;
526     const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
527     const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
528     const index_t N2=(neighbour/(m_NX*m_NY)==0 ? m_N2 : m_N2-1);
529     #pragma omp parallel for
530     for (dim_t i1=bottom; i1<m_N1; i1++) {
531     m_nodeId[i1*m_N0]=m_nodeDistribution[neighbour]
532     + N0*N1*(N2-1)+N0-1+(i1-bottom)*N0;
533     }
534     }
535     // case 7: bottom-left-front corner node owned by corresponding rank
536     if (left>0 && bottom>0 && front>0) {
537     const int neighbour=m_mpiInfo->rank-m_NX*(m_NY+1)-1;
538     const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
539     const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
540     const index_t N2=(neighbour/(m_NX*m_NY) == 0 ? m_N2 : m_N2-1);
541     m_nodeId[0]=m_nodeDistribution[neighbour]+N0*N1*N2-1;
542     }
543 caltinay 3697
544 caltinay 3698 // the rest of the id's are contiguous
545     const index_t firstId=m_nodeDistribution[m_mpiInfo->rank];
546     #pragma omp parallel for
547     for (dim_t i2=front; i2<m_N2; i2++) {
548     for (dim_t i1=bottom; i1<m_N1; i1++) {
549     for (dim_t i0=left; i0<m_N0; i0++) {
550     m_nodeId[i0+i1*m_N0+i2*m_N0*m_N1] = firstId+i0-left
551     +(i1-bottom)*(m_N0-left)
552     +(i2-front)*(m_N0-left)*(m_N1-bottom);
553     }
554     }
555     }
556    
557 caltinay 3697 // elements
558     m_elementId.resize(getNumElements());
559     #pragma omp parallel for
560     for (dim_t k=0; k<getNumElements(); k++) {
561     m_elementId[k]=k;
562     }
563    
564     // face elements
565     m_faceId.resize(getNumFaceElements());
566     #pragma omp parallel for
567     for (dim_t k=0; k<getNumFaceElements(); k++) {
568     m_faceId[k]=k;
569     }
570 caltinay 3691 }
571    
572     } // end of namespace ripley
573    

  ViewVC Help
Powered by ViewVC 1.1.26