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

Diff of /branches/diaplayground/ripley/src/Brick.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3697 by caltinay, Mon Nov 28 04:52:00 2011 UTC revision 3698 by caltinay, Tue Nov 29 00:47:29 2011 UTC
# Line 147  void Brick::dump(const string& fileName) Line 147  void Brick::dump(const string& fileName)
147      boost::scoped_ptr<double> y(new double[m_N1]);      boost::scoped_ptr<double> y(new double[m_N1]);
148      boost::scoped_ptr<double> z(new double[m_N2]);      boost::scoped_ptr<double> z(new double[m_N2]);
149      double* coords[3] = { x.get(), y.get(), z.get() };      double* coords[3] = { x.get(), y.get(), z.get() };
150        pair<double,double> xdx = getFirstCoordAndSpacing(0);
151        pair<double,double> ydy = getFirstCoordAndSpacing(1);
152        pair<double,double> zdz = getFirstCoordAndSpacing(2);
153  #pragma omp parallel  #pragma omp parallel
154      {      {
155  #pragma omp for  #pragma omp for
156          for (dim_t i0 = 0; i0 < m_N0; i0++) {          for (dim_t i0 = 0; i0 < m_N0; i0++) {
157              coords[0][i0]=(m_l0*(i0+m_offset0))/m_gNE0;              coords[0][i0]=xdx.first+i0*xdx.second;
158          }          }
159  #pragma omp for  #pragma omp for
160          for (dim_t i1 = 0; i1 < m_N1; i1++) {          for (dim_t i1 = 0; i1 < m_N1; i1++) {
161              coords[1][i1]=(m_l1*(i1+m_offset1))/m_gNE1;              coords[1][i1]=ydy.first+i1*ydy.second;
162          }          }
163  #pragma omp for  #pragma omp for
164          for (dim_t i2 = 0; i2 < m_N2; i2++) {          for (dim_t i2 = 0; i2 < m_N2; i2++) {
165              coords[2][i2]=(m_l2*(i2+m_offset2))/m_gNE2;              coords[2][i2]=zdz.first+i2*zdz.second;
166          }          }
167      }      }
168      int dims[] = { m_N0, m_N1, m_N2 };      IndexVector dims = getNumNodesPerDim();
169      DBPutQuadmesh(dbfile, "mesh", NULL, coords, dims, 3, DB_DOUBLE,      DBPutQuadmesh(dbfile, "mesh", NULL, coords, &dims[0], 3, DB_DOUBLE,
170              DB_COLLINEAR, NULL);              DB_COLLINEAR, NULL);
171    
172      DBPutQuadvar1(dbfile, "nodeId", "mesh", (void*)&m_nodeId[0], dims, 3, NULL,      DBPutQuadvar1(dbfile, "nodeId", "mesh", (void*)&m_nodeId[0], &dims[0], 3,
173              0, DB_INT, DB_NODECENT, NULL);              NULL, 0, DB_INT, DB_NODECENT, NULL);
174    
175        // 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      if (m_mpiInfo->rank == 0) {      if (m_mpiInfo->rank == 0) {
182          vector<string> tempstrings;          vector<string> tempstrings;
183          vector<char*> names;          vector<char*> names;
# Line 193  void Brick::dump(const string& fileName) Line 202  void Brick::dump(const string& fileName)
202          types.assign(m_mpiInfo->size, DB_QUADVAR);          types.assign(m_mpiInfo->size, DB_QUADVAR);
203          DBPutMultivar(dbfile, "nodeId", m_mpiInfo->size, &names[0],          DBPutMultivar(dbfile, "nodeId", m_mpiInfo->size, &names[0],
204                 &types[0], NULL);                 &types[0], NULL);
205            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      }      }
216    
217      if (m_mpiInfo->size > 1) {      if (m_mpiInfo->size > 1) {
# Line 232  bool Brick::ownSample(int fsCode, index_ Line 251  bool Brick::ownSample(int fsCode, index_
251  {  {
252  #ifdef ESYS_MPI  #ifdef ESYS_MPI
253      if (fsCode == Nodes) {      if (fsCode == Nodes) {
254          const index_t myFirst=getNumNodes()*m_mpiInfo->rank;          const index_t myFirst=m_nodeDistribution[m_mpiInfo->rank];
255          const index_t myLast=getNumNodes()*(m_mpiInfo->rank+1)-1;          const index_t myLast=m_nodeDistribution[m_mpiInfo->rank+1]-1;
256          return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);          return (m_nodeId[id]>=myFirst && m_nodeId[id]<=myLast);
257      } else      } else
258          throw RipleyException("ownSample() only implemented for Nodes");          throw RipleyException("ownSample() only implemented for Nodes");
# Line 271  void Brick::Print_Mesh_Info(const bool f Line 290  void Brick::Print_Mesh_Info(const bool f
290          cout << "     Id  Coordinates" << endl;          cout << "     Id  Coordinates" << endl;
291          cout.precision(15);          cout.precision(15);
292          cout.setf(ios::scientific, ios::floatfield);          cout.setf(ios::scientific, ios::floatfield);
293            pair<double,double> xdx = getFirstCoordAndSpacing(0);
294            pair<double,double> ydy = getFirstCoordAndSpacing(1);
295            pair<double,double> zdz = getFirstCoordAndSpacing(2);
296          for (index_t i=0; i < getNumNodes(); i++) {          for (index_t i=0; i < getNumNodes(); i++) {
297              cout << "  " << setw(5) << m_nodeId[i]              cout << "  " << setw(5) << m_nodeId[i]
298                  << "  " << (m_l0*(i%m_N0+m_offset0))/m_gNE0                  << "  " << xdx.first+(i%m_N0)*xdx.second
299                  << "  " << (m_l1*(i%(m_N0*m_N1)/m_N0+m_offset1))/m_gNE1                  << "  " << ydy.first+(i%(m_N0*m_N1)/m_N0)*ydy.second
300                  << "  " << (m_l2*(i/(m_N0*m_N1)+m_offset2))/m_gNE2 << endl;                  << "  " << zdz.first+(i/(m_N0*m_N1))*zdz.second << endl;
301          }          }
302      }      }
303  }  }
304    
305    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  //protected  //protected
361  dim_t Brick::getNumFaceElements() const  dim_t Brick::getNumFaceElements() const
362  {  {
# Line 316  void Brick::assembleCoordinates(escript: Line 393  void Brick::assembleCoordinates(escript:
393      if (!numSamplesEqual(&x, 1, getNumNodes()))      if (!numSamplesEqual(&x, 1, getNumNodes()))
394          throw RipleyException("setToX: Illegal number of samples in Data object");          throw RipleyException("setToX: Illegal number of samples in Data object");
395    
396        pair<double,double> xdx = getFirstCoordAndSpacing(0);
397        pair<double,double> ydy = getFirstCoordAndSpacing(1);
398        pair<double,double> zdz = getFirstCoordAndSpacing(2);
399      arg.requireWrite();      arg.requireWrite();
400  #pragma omp parallel for  #pragma omp parallel for
401      for (dim_t i2 = 0; i2 < m_N2; i2++) {      for (dim_t i2 = 0; i2 < m_N2; i2++) {
402          for (dim_t i1 = 0; i1 < m_N1; i1++) {          for (dim_t i1 = 0; i1 < m_N1; i1++) {
403              for (dim_t i0 = 0; i0 < m_N0; i0++) {              for (dim_t i0 = 0; i0 < m_N0; i0++) {
404                  double* point = arg.getSampleDataRW(i0+m_N0*i1+m_N0*m_N1*i2);                  double* point = arg.getSampleDataRW(i0+m_N0*i1+m_N0*m_N1*i2);
405                  point[0] = (m_l0*(i0+m_offset0))/m_gNE0;                  point[0] = xdx.first+i0*xdx.second;
406                  point[1] = (m_l1*(i1+m_offset1))/m_gNE1;                  point[1] = ydy.first+i1*ydy.second;
407                  point[2] = (m_l2*(i2+m_offset2))/m_gNE2;                  point[2] = zdz.first+i2*zdz.second;
408              }              }
409          }          }
410      }      }
# Line 336  void Brick::populateSampleIds() Line 416  void Brick::populateSampleIds()
416      // identifiers are ordered from left to right, bottom to top, front to back      // 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      // 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      // below / to the left / to the front of the current rank
419      const index_t firstId = getNumNodes()*m_mpiInfo->rank;  
420      const index_t diff0 = m_N0*(m_N1*m_N2-1)+1;      // build node distribution vector first.
421      const index_t diff1 = m_N0*m_N1*(m_N2*m_NX-1)+m_N0;      // m_nodeDistribution[i] is the first node id on rank i, that is
422      const index_t diff2 = m_N0*m_N1*m_N2*m_NX*m_NY-m_N0*m_N1*(m_N2-1);      // 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      m_nodeId.resize(getNumNodes());      m_nodeId.resize(getNumNodes());
451    
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    #pragma omp parallel for
464            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        }
472        // 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    
544        // the rest of the id's are contiguous
545        const index_t firstId=m_nodeDistribution[m_mpiInfo->rank];
546  #pragma omp parallel for  #pragma omp parallel for
547      for (dim_t k=0; k<getNumNodes(); k++) {      for (dim_t i2=front; i2<m_N2; i2++) {
548          index_t id = firstId+k;          for (dim_t i1=bottom; i1<m_N1; i1++) {
549          if (m_offset0 > 0 && k%m_N0==0)              for (dim_t i0=left; i0<m_N0; i0++) {
550              id -= diff0;                  m_nodeId[i0+i1*m_N0+i2*m_N0*m_N1] = firstId+i0-left
551          if (m_offset1 > 0 && k%(m_N0*m_N1)<m_N0)                      +(i1-bottom)*(m_N0-left)
552              id -= diff1;                      +(i2-front)*(m_N0-left)*(m_N1-bottom);
553          if (m_offset2 > 0 && k/(m_N0*m_N1)==0)              }
554              id -= diff2;          }
         m_nodeId[k]=id;  
555      }      }
556    
557      // elements      // elements

Legend:
Removed from v.3697  
changed lines
  Added in v.3698

  ViewVC Help
Powered by ViewVC 1.1.26