/[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 3691 by caltinay, Wed Nov 23 23:07:37 2011 UTC revision 3702 by caltinay, Fri Dec 2 06:12:32 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, 0,      DBPutQuadvar1(dbfile, "nodeId", "mesh", (void*)&m_nodeId[0], &dims[0], 3,
173              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 211  void Brick::dump(const string& fileName) Line 230  void Brick::dump(const string& fileName)
230    
231  const int* Brick::borrowSampleReferenceIDs(int fsType) const  const int* Brick::borrowSampleReferenceIDs(int fsType) const
232  {  {
233      if (fsType == Nodes)      switch (fsType) {
234          return &m_nodeId[0];          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    
244      throw RipleyException("borrowSampleReferenceIDs() only implemented for Nodes");      stringstream msg;
245        msg << "borrowSampleReferenceIDs() not implemented for function space type "
246            << fsType;
247        throw RipleyException(msg.str());
248  }  }
249    
250  bool Brick::ownSample(int fsCode, index_t id) const  bool Brick::ownSample(int fsCode, index_t id) const
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 231  bool Brick::ownSample(int fsCode, index_ Line 261  bool Brick::ownSample(int fsCode, index_
261  #endif  #endif
262  }  }
263    
 void Brick::interpolateOnDomain(escript::Data& target,  
                                 const escript::Data& in) const  
 {  
     const Brick& inDomain=dynamic_cast<const Brick&>(*(in.getFunctionSpace().getDomain()));  
     const Brick& targetDomain=dynamic_cast<const Brick&>(*(target.getFunctionSpace().getDomain()));  
     if (inDomain != *this)  
         throw RipleyException("Illegal domain of interpolant");  
     if (targetDomain != *this)  
         throw RipleyException("Illegal domain of interpolation target");  
   
     throw RipleyException("interpolateOnDomain() not implemented");  
 }  
   
264  Paso_SystemMatrixPattern* Brick::getPattern(bool reducedRowOrder,  Paso_SystemMatrixPattern* Brick::getPattern(bool reducedRowOrder,
265                                              bool reducedColOrder) const                                              bool reducedColOrder) const
266  {  {
# Line 260  void Brick::Print_Mesh_Info(const bool f Line 277  void Brick::Print_Mesh_Info(const bool f
277          cout << "     Id  Coordinates" << endl;          cout << "     Id  Coordinates" << endl;
278          cout.precision(15);          cout.precision(15);
279          cout.setf(ios::scientific, ios::floatfield);          cout.setf(ios::scientific, ios::floatfield);
280            pair<double,double> xdx = getFirstCoordAndSpacing(0);
281            pair<double,double> ydy = getFirstCoordAndSpacing(1);
282            pair<double,double> zdz = getFirstCoordAndSpacing(2);
283          for (index_t i=0; i < getNumNodes(); i++) {          for (index_t i=0; i < getNumNodes(); i++) {
284              cout << "  " << setw(5) << m_nodeId[i]              cout << "  " << setw(5) << m_nodeId[i]
285                  << "  " << (m_l0*(i%m_N0+m_offset0))/m_gNE0                  << "  " << xdx.first+(i%m_N0)*xdx.second
286                  << "  " << (m_l1*(i%(m_N0*m_N1)/m_N0+m_offset1))/m_gNE1                  << "  " << ydy.first+(i%(m_N0*m_N1)/m_N0)*ydy.second
287                  << "  " << (m_l2*(i/(m_N0*m_N1)+m_offset2))/m_gNE2 << endl;                  << "  " << zdz.first+(i/(m_N0*m_N1))*zdz.second << endl;
288          }          }
289      }      }
290  }  }
291    
292    IndexVector Brick::getNumNodesPerDim() const
293    {
294        IndexVector ret;
295        ret.push_back(m_N0);
296        ret.push_back(m_N1);
297        ret.push_back(m_N2);
298        return ret;
299    }
300    
301    IndexVector Brick::getNumElementsPerDim() const
302    {
303        IndexVector ret;
304        ret.push_back(m_NE0);
305        ret.push_back(m_NE1);
306        ret.push_back(m_NE2);
307        return ret;
308    }
309    
310    IndexVector Brick::getNumFacesPerBoundary() const
311    {
312        IndexVector ret(6, 0);
313        //left
314        if (m_offset0==0)
315            ret[0]=m_NE1*m_NE2;
316        //right
317        if (m_mpiInfo->rank%m_NX==m_NX-1)
318            ret[1]=m_NE1*m_NE2;
319        //bottom
320        if (m_offset1==0)
321            ret[2]=m_NE0*m_NE2;
322        //top
323        if (m_mpiInfo->rank%(m_NX*m_NY)/m_NX==m_NY-1)
324            ret[3]=m_NE0*m_NE2;
325        //front
326        if (m_offset2==0)
327            ret[4]=m_NE0*m_NE1;
328        //back
329        if (m_mpiInfo->rank/(m_NX*m_NY)==m_NZ-1)
330            ret[5]=m_NE0*m_NE1;
331        return ret;
332    }
333    
334    pair<double,double> Brick::getFirstCoordAndSpacing(dim_t dim) const
335    {
336        if (dim==0)
337            return pair<double,double>((m_l0*m_offset0)/m_gNE0, m_l0/m_gNE0);
338        else if (dim==1)
339            return pair<double,double>((m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);
340        else if (dim==2)
341            return pair<double,double>((m_l2*m_offset2)/m_gNE2, m_l2/m_gNE2);
342    
343        throw RipleyException("getFirstCoordAndSpacing(): invalid argument");
344    }
345    
346    
347  //protected  //protected
348  dim_t Brick::getNumFaceElements() const  dim_t Brick::getNumFaceElements() const
349  {  {
# Line 305  void Brick::assembleCoordinates(escript: Line 380  void Brick::assembleCoordinates(escript:
380      if (!numSamplesEqual(&x, 1, getNumNodes()))      if (!numSamplesEqual(&x, 1, getNumNodes()))
381          throw RipleyException("setToX: Illegal number of samples in Data object");          throw RipleyException("setToX: Illegal number of samples in Data object");
382    
383        pair<double,double> xdx = getFirstCoordAndSpacing(0);
384        pair<double,double> ydy = getFirstCoordAndSpacing(1);
385        pair<double,double> zdz = getFirstCoordAndSpacing(2);
386      arg.requireWrite();      arg.requireWrite();
387  #pragma omp parallel for  #pragma omp parallel for
388      for (dim_t i2 = 0; i2 < m_N2; i2++) {      for (dim_t i2 = 0; i2 < m_N2; i2++) {
389          for (dim_t i1 = 0; i1 < m_N1; i1++) {          for (dim_t i1 = 0; i1 < m_N1; i1++) {
390              for (dim_t i0 = 0; i0 < m_N0; i0++) {              for (dim_t i0 = 0; i0 < m_N0; i0++) {
391                  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);
392                  point[0] = (m_l0*(i0+m_offset0))/m_gNE0;                  point[0] = xdx.first+i0*xdx.second;
393                  point[1] = (m_l1*(i1+m_offset1))/m_gNE1;                  point[1] = ydy.first+i1*ydy.second;
394                  point[2] = (m_l2*(i2+m_offset2))/m_gNE2;                  point[2] = zdz.first+i2*zdz.second;
395              }              }
396          }          }
397      }      }
# Line 322  void Brick::assembleCoordinates(escript: Line 400  void Brick::assembleCoordinates(escript:
400  //private  //private
401  void Brick::populateSampleIds()  void Brick::populateSampleIds()
402  {  {
403      const index_t firstId = getNumNodes()*m_mpiInfo->rank;      // identifiers are ordered from left to right, bottom to top, front to back
404      const index_t diff0 = m_N0*(m_N1*m_N2-1)+1;      // on each rank, except for the shared nodes which are owned by the rank
405      const index_t diff1 = m_N0*m_N1*(m_N2*m_NX-1)+m_N0;      // below / to the left / to the front of the current rank
406      const index_t diff2 = m_N0*m_N1*m_N2*m_NX*m_NY-m_N0*m_N1*(m_N2-1);  
407        // build node distribution vector first.
408        // m_nodeDistribution[i] is the first node id on rank i, that is
409        // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes
410        m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
411        m_nodeDistribution[1]=getNumNodes();
412        for (dim_t k=1; k<m_mpiInfo->size-1; k++) {
413            const index_t x = k%m_NX;
414            const index_t y = k%(m_NX*m_NY)/m_NX;
415            const index_t z = k/(m_NX*m_NY);
416            index_t numNodes=getNumNodes();
417            if (x>0)
418                numNodes-=m_N1*m_N2;
419            if (y>0)
420                numNodes-=m_N0*m_N2;
421            if (z>0)
422                numNodes-=m_N0*m_N1;
423            // if an edge was subtracted twice add it back
424            if (x>0 && y>0)
425                numNodes+=m_N2;
426            if (x>0 && z>0)
427                numNodes+=m_N1;
428            if (y>0 && z>0)
429                numNodes+=m_N0;
430            // the corner node was removed 3x and added back 3x, so subtract it
431            if (x>0 && y>0 && z>0)
432                numNodes--;
433            m_nodeDistribution[k+1]=m_nodeDistribution[k]+numNodes;
434        }
435        m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
436    
437      m_nodeId.resize(getNumNodes());      m_nodeId.resize(getNumNodes());
438    
439        // the bottom, left and front planes are not owned by this rank so the
440        // identifiers need to be computed accordingly
441        const index_t left = (m_offset0==0 ? 0 : 1);
442        const index_t bottom = (m_offset1==0 ? 0 : 1);
443        const index_t front = (m_offset2==0 ? 0 : 1);
444    
445        // case 1: all nodes on left plane are owned by rank on the left
446        if (left>0) {
447            const int neighbour=m_mpiInfo->rank-1;
448            const index_t leftN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
449            const index_t leftN1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
450    #pragma omp parallel for
451            for (dim_t i2=front; i2<m_N2; i2++) {
452                for (dim_t i1=bottom; i1<m_N1; i1++) {
453                    m_nodeId[i1*m_N0+i2*m_N0*m_N1]=m_nodeDistribution[neighbour]
454                        + (i1-bottom+1)*leftN0
455                        + (i2-front)*leftN0*leftN1 - 1;
456                }
457            }
458        }
459        // case 2: all nodes on bottom plane are owned by rank below
460        if (bottom>0) {
461            const int neighbour=m_mpiInfo->rank-m_NX;
462            const index_t bottomN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
463            const index_t bottomN1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
464    #pragma omp parallel for
465            for (dim_t i2=front; i2<m_N2; i2++) {
466                for (dim_t i0=left; i0<m_N0; i0++) {
467                    m_nodeId[i0+i2*m_N0*m_N1]=m_nodeDistribution[neighbour]
468                        + bottomN0*(bottomN1-1)
469                        + (i2-front)*bottomN0*bottomN1 + i0-left;
470                }
471            }
472        }
473        // case 3: all nodes on front plane are owned by rank in front
474        if (front>0) {
475            const int neighbour=m_mpiInfo->rank-m_NX*m_NY;
476            const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
477            const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
478            const index_t N2=(neighbour/(m_NX*m_NY)==0 ? m_N2 : m_N2-1);
479    #pragma omp parallel for
480            for (dim_t i1=bottom; i1<m_N1; i1++) {
481                for (dim_t i0=left; i0<m_N0; i0++) {
482                    m_nodeId[i0+i1*m_N0]=m_nodeDistribution[neighbour]
483                        + N0*N1*(N2-1)+(i1-bottom)*N0 + i0-left;
484                }
485            }
486        }
487        // case 4: nodes on front bottom edge are owned by the corresponding rank
488        if (front>0 && bottom>0) {
489            const int neighbour=m_mpiInfo->rank-m_NX*(m_NY+1);
490            const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
491            const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
492            const index_t N2=(neighbour/(m_NX*m_NY)==0 ? m_N2 : m_N2-1);
493    #pragma omp parallel for
494            for (dim_t i0=left; i0<m_N0; i0++) {
495                m_nodeId[i0]=m_nodeDistribution[neighbour]
496                    + N0*N1*(N2-1)+(N1-1)*N0 + i0-left;
497            }
498        }
499        // case 5: nodes on left bottom edge are owned by the corresponding rank
500        if (left>0 && bottom>0) {
501            const int neighbour=m_mpiInfo->rank-m_NX-1;
502            const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
503            const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
504    #pragma omp parallel for
505            for (dim_t i2=front; i2<m_N2; i2++) {
506                m_nodeId[i2*m_N0*m_N1]=m_nodeDistribution[neighbour]
507                    + (1+i2-front)*N0*N1-1;
508            }
509        }
510        // case 6: nodes on left front edge are owned by the corresponding rank
511        if (left>0 && front>0) {
512            const int neighbour=m_mpiInfo->rank-m_NX*m_NY-1;
513            const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
514            const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
515            const index_t N2=(neighbour/(m_NX*m_NY)==0 ? m_N2 : m_N2-1);
516    #pragma omp parallel for
517            for (dim_t i1=bottom; i1<m_N1; i1++) {
518                m_nodeId[i1*m_N0]=m_nodeDistribution[neighbour]
519                    + N0*N1*(N2-1)+N0-1+(i1-bottom)*N0;
520            }
521        }
522        // case 7: bottom-left-front corner node owned by corresponding rank
523        if (left>0 && bottom>0 && front>0) {
524            const int neighbour=m_mpiInfo->rank-m_NX*(m_NY+1)-1;
525            const index_t N0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
526            const index_t N1=(neighbour%(m_NX*m_NY)/m_NX==0 ? m_N1 : m_N1-1);
527            const index_t N2=(neighbour/(m_NX*m_NY) == 0 ? m_N2 : m_N2-1);
528            m_nodeId[0]=m_nodeDistribution[neighbour]+N0*N1*N2-1;
529        }
530    
531        // the rest of the id's are contiguous
532        const index_t firstId=m_nodeDistribution[m_mpiInfo->rank];
533    #pragma omp parallel for
534        for (dim_t i2=front; i2<m_N2; i2++) {
535            for (dim_t i1=bottom; i1<m_N1; i1++) {
536                for (dim_t i0=left; i0<m_N0; i0++) {
537                    m_nodeId[i0+i1*m_N0+i2*m_N0*m_N1] = firstId+i0-left
538                        +(i1-bottom)*(m_N0-left)
539                        +(i2-front)*(m_N0-left)*(m_N1-bottom);
540                }
541            }
542        }
543    
544        // elements
545        m_elementId.resize(getNumElements());
546    #pragma omp parallel for
547        for (dim_t k=0; k<getNumElements(); k++) {
548            m_elementId[k]=k;
549        }
550    
551        // face elements
552        m_faceId.resize(getNumFaceElements());
553  #pragma omp parallel for  #pragma omp parallel for
554      for (dim_t k=0; k<getNumNodes(); k++) {      for (dim_t k=0; k<getNumFaceElements(); k++) {
555          index_t id = firstId+k;          m_faceId[k]=k;
         if (m_offset0 > 0 && k%m_N0==0)  
             id -= diff0;  
         if (m_offset1 > 0 && k%(m_N0*m_N1)<m_N0)  
             id -= diff1;  
         if (m_offset2 > 0 && k/(m_N0*m_N1)==0)  
             id -= diff2;  
         m_nodeId[k]=id;  
556      }      }
557  }  }
558    

Legend:
Removed from v.3691  
changed lines
  Added in v.3702

  ViewVC Help
Powered by ViewVC 1.1.26