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

Diff of /trunk/ripley/src/Rectangle.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 3697 by caltinay, Mon Nov 28 04:52:00 2011 UTC
# Line 138  void Rectangle::dump(const string& fileN Line 138  void Rectangle::dump(const string& fileN
138      boost::scoped_ptr<double> x(new double[m_N0]);      boost::scoped_ptr<double> x(new double[m_N0]);
139      boost::scoped_ptr<double> y(new double[m_N1]);      boost::scoped_ptr<double> y(new double[m_N1]);
140      double* coords[2] = { x.get(), y.get() };      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  #pragma omp parallel
144      {      {
145  #pragma omp for  #pragma omp for
146          for (dim_t i0 = 0; i0 < m_N0; i0++) {          for (dim_t i0 = 0; i0 < m_N0; i0++) {
147              coords[0][i0]=(m_l0*(i0+m_offset0))/m_gNE0;              coords[0][i0]=xdx.first+i0*xdx.second;
148          }          }
149  #pragma omp for  #pragma omp for
150          for (dim_t i1 = 0; i1 < m_N1; i1++) {          for (dim_t i1 = 0; i1 < m_N1; i1++) {
151              coords[1][i1]=(m_l1*(i1+m_offset1))/m_gNE1;              coords[1][i1]=ydy.first+i1*ydy.second;
152          }          }
153      }      }
154      int dims[] = { m_N0, m_N1 };      IndexVector dims = getNumNodesPerDim();
155      DBPutQuadmesh(dbfile, "mesh", NULL, coords, dims, 2, DB_DOUBLE,  
156        // write mesh
157        DBPutQuadmesh(dbfile, "mesh", NULL, coords, &dims[0], 2, DB_DOUBLE,
158              DB_COLLINEAR, NULL);              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) {      if (m_mpiInfo->rank == 0) {
171          vector<string> tempstrings;          vector<string> tempstrings;
172          vector<char*> meshnames;          vector<char*> names;
173          for (dim_t i=0; i<m_mpiInfo->size; i++) {          for (dim_t i=0; i<m_mpiInfo->size; i++) {
174              stringstream path;              stringstream path;
175              path << "/block" << setw(4) << setfill('0') << right << i << "/mesh";              path << "/block" << setw(4) << setfill('0') << right << i << "/mesh";
176              tempstrings.push_back(path.str());              tempstrings.push_back(path.str());
177              meshnames.push_back((char*)tempstrings.back().c_str());              names.push_back((char*)tempstrings.back().c_str());
178          }          }
179          vector<int> meshtypes(m_mpiInfo->size, DB_QUAD_RECT);          vector<int> types(m_mpiInfo->size, DB_QUAD_RECT);
180          DBSetDir(dbfile, "/");          DBSetDir(dbfile, "/");
181          DBPutMultimesh(dbfile, "multimesh", m_mpiInfo->size, &meshnames[0],          DBPutMultimesh(dbfile, "multimesh", m_mpiInfo->size, &names[0],
182                 &meshtypes[0], NULL);                 &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) {      if (m_mpiInfo->size > 1) {
# Line 182  void Rectangle::dump(const string& fileN Line 217  void Rectangle::dump(const string& fileN
217  #endif  #endif
218  }  }
219    
220  const int* Rectangle::borrowSampleReferenceIDs(int functionSpaceType) const  const int* Rectangle::borrowSampleReferenceIDs(int fsType) const
221  {  {
222      switch (functionSpaceType) {      switch (fsType) {
223          case Nodes:          case Nodes:
224              return &m_nodeId[0];              return &m_nodeId[0];
225          case Elements:          case Elements:
# Line 197  const int* Rectangle::borrowSampleRefere Line 232  const int* Rectangle::borrowSampleRefere
232    
233      stringstream msg;      stringstream msg;
234      msg << "borrowSampleReferenceIDs() not implemented for function space type "      msg << "borrowSampleReferenceIDs() not implemented for function space type "
235          << functionSpaceType;          << fsType;
236      throw RipleyException(msg.str());      throw RipleyException(msg.str());
237  }  }
238    
# Line 225  void Rectangle::interpolateOnDomain(escr Line 260  void Rectangle::interpolateOnDomain(escr
260      if (targetDomain != *this)      if (targetDomain != *this)
261          throw RipleyException("Illegal domain of interpolation target");          throw RipleyException("Illegal domain of interpolation target");
262    
263      throw RipleyException("interpolateOnDomain() not implemented");      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 Elements:
273                    {
274                        const double tmp0_2 = 0.62200846792814621559;
275                        const double tmp0_1 = 0.044658198738520451079;
276                        const double tmp0_0 = 0.16666666666666666667;
277                        const dim_t numComp = in.getDataPointSize();
278                        escript::Data* inn=const_cast<escript::Data*>(&in);
279    #pragma omp parallel for
280                        for (index_t k1=0; k1 < m_NE1; ++k1) {
281                            for (index_t k0=0; k0 < m_NE0; ++k0) {
282                                const register double* f_10 = inn->getSampleDataRO(INDEX2(k0+1,k1, m_N0));
283                                const register double* f_11 = inn->getSampleDataRO(INDEX2(k0+1,k1+1, m_N0));
284                                const register double* f_01 = inn->getSampleDataRO(INDEX2(k0,k1+1, m_N0));
285                                const register double* f_00 = inn->getSampleDataRO(INDEX2(k0,k1, m_N0));
286                                double* o = target.getSampleDataRW(INDEX2(k0,k1,m_NE0));
287                                for (index_t i=0; i < numComp; ++i) {
288                                    o[INDEX2(i,numComp,0)] = f_00[i]*tmp0_2 + f_11[i]*tmp0_1 + tmp0_0*(f_01[i] + f_10[i]);
289                                    o[INDEX2(i,numComp,1)] = f_01[i]*tmp0_1 + f_10[i]*tmp0_2 + tmp0_0*(f_00[i] + f_11[i]);
290                                    o[INDEX2(i,numComp,2)] = f_01[i]*tmp0_2 + f_10[i]*tmp0_1 + tmp0_0*(f_00[i] + f_11[i]);
291                                    o[INDEX2(i,numComp,3)] = f_00[i]*tmp0_1 + f_11[i]*tmp0_2 + tmp0_0*(f_01[i] + f_10[i]);
292                                } // close component loop i
293                            } // close k0 loop
294                        } // close k1 loop
295                        break;
296                    }
297    
298                    case DegreesOfFreedom:
299                        target=in;
300                        break;
301                        
302                    default:
303                        throw RipleyException(msg.str());
304                }
305            default:
306                throw RipleyException(msg.str());
307        }
308  }  }
309    
310  Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,  Paso_SystemMatrixPattern* Rectangle::getPattern(bool reducedRowOrder,
# Line 243  Paso_SystemMatrixPattern* Rectangle::get Line 322  Paso_SystemMatrixPattern* Rectangle::get
322              &dist[0], 1, 0);              &dist[0], 1, 0);
323    
324      // connectors      // connectors
325      dim_t numNeighbours = ...;      dim_t numNeighbours = 0;
326      RankVector neighbour(numNeighbours);      RankVector neighbour(numNeighbours);
327      IndexVector offsetInShared(numNeighbours+1);      IndexVector offsetInShared(numNeighbours+1);
328      IndexVector shared(offsetInShared[numNeighbours]);      IndexVector shared(offsetInShared[numNeighbours]);
329    
330      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(      Paso_SharedComponents *snd_shcomp = Paso_SharedComponents_alloc(
331              getNumNodes(), numNeighbours, neighbour, shared, offsetInShared,              getNumNodes(), numNeighbours, &neighbour[0], &shared[0],
332              1, 0, m_mpiInfo);              &offsetInShared[0], 1, 0, m_mpiInfo);
333      Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(      Paso_SharedComponents *rcv_shcomp = Paso_SharedComponents_alloc(
334              getNumNodes(), numNeighbours, neighbour, shared, offsetInShared,              getNumNodes(), numNeighbours, &neighbour[0], &shared[0],
335              1, 0, m_mpiInfo);              &offsetInShared[0], 1, 0, m_mpiInfo);
336      Paso_Connector* connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);      Paso_Connector* connector = Paso_Connector_alloc(snd_shcomp, rcv_shcomp);
337    
338      // create patterns      // create patterns
339        dim_t M=0, N=0;
340        int* ptr=NULL;
341        index_t* index=NULL;
342      Paso_Pattern *mainPattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,      Paso_Pattern *mainPattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,
343              M, N, ptr, index);              M, N, ptr, index);
344      Paso_Pattern *colCouplePattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,      Paso_Pattern *colCouplePattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,
345              M, N, ...);              M, N, ptr, index);
346      Paso_Pattern *rowCouplePattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,      Paso_Pattern *rowCouplePattern = Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,
347              ...);              M, N, ptr, index);
348    
349      Paso_SystemMatrixPattern* pattern = Paso_SystemMatrixPattern_alloc(      Paso_SystemMatrixPattern* pattern = Paso_SystemMatrixPattern_alloc(
350              MATRIX_FORMAT_DEFAULT, distribution, distribution,              MATRIX_FORMAT_DEFAULT, distribution, distribution,
# Line 274  Paso_SystemMatrixPattern* Rectangle::get Line 356  Paso_SystemMatrixPattern* Rectangle::get
356      Paso_Distribution_free(distribution);      Paso_Distribution_free(distribution);
357      Paso_SharedComponents_free(snd_shcomp);      Paso_SharedComponents_free(snd_shcomp);
358      Paso_SharedComponents_free(rcv_shcomp);      Paso_SharedComponents_free(rcv_shcomp);
359        return pattern;
360  */  */
361      throw RipleyException("getPattern() not implemented");      throw RipleyException("getPattern() not implemented");
362  }  }
# Line 285  void Rectangle::Print_Mesh_Info(const bo Line 368  void Rectangle::Print_Mesh_Info(const bo
368          cout << "     Id  Coordinates" << endl;          cout << "     Id  Coordinates" << endl;
369          cout.precision(15);          cout.precision(15);
370          cout.setf(ios::scientific, ios::floatfield);          cout.setf(ios::scientific, ios::floatfield);
371            pair<double,double> xdx = getFirstCoordAndSpacing(0);
372            pair<double,double> ydy = getFirstCoordAndSpacing(1);
373          for (index_t i=0; i < getNumNodes(); i++) {          for (index_t i=0; i < getNumNodes(); i++) {
374              cout << "  " << setw(5) << m_nodeId[i]              cout << "  " << setw(5) << m_nodeId[i]
375                  << "  " << (m_l0*(i%m_N0+m_offset0))/m_gNE0                  << "  " << xdx.first+(i%m_N0)*xdx.second
376                  << "  " << (m_l1*(i/m_N0+m_offset1))/m_gNE1 << endl;                  << "  " << ydy.first+(i/m_N0)*ydy.second << endl;
377          }          }
378      }      }
379  }  }
380    
381    IndexVector Rectangle::getNumNodesPerDim() const
382    {
383        IndexVector ret;
384        ret.push_back(m_N0);
385        ret.push_back(m_N1);
386        return ret;
387    }
388    
389    IndexVector Rectangle::getNumElementsPerDim() const
390    {
391        IndexVector ret;
392        ret.push_back(m_NE0);
393        ret.push_back(m_NE1);
394        return ret;
395    }
396    
397    IndexVector Rectangle::getNumFacesPerBoundary() const
398    {
399        IndexVector ret(4, 0);
400        //left
401        if (m_offset0==0)
402            ret[0]=m_NE1;
403        //right
404        if (m_mpiInfo->rank%m_NX==m_NX-1)
405            ret[1]=m_NE1;
406        //bottom
407        if (m_offset1==0)
408            ret[2]=m_NE0;
409        //top
410        if (m_mpiInfo->rank/m_NX==m_NY-1)
411            ret[3]=m_NE0;
412        return ret;
413    }
414    
415    pair<double,double> Rectangle::getFirstCoordAndSpacing(dim_t dim) const
416    {
417        if (dim==0) {
418            return pair<double,double>((m_l0*m_offset0)/m_gNE0, m_l0/m_gNE0);
419        } else if (dim==1) {
420            return pair<double,double>((m_l1*m_offset1)/m_gNE1, m_l1/m_gNE1);
421        }
422        throw RipleyException("getFirstCoordAndSpacing(): invalid argument");
423    }
424    
425  //protected  //protected
426  dim_t Rectangle::getNumFaceElements() const  dim_t Rectangle::getNumFaceElements() const
427  {  {
# Line 323  void Rectangle::assembleCoordinates(escr Line 452  void Rectangle::assembleCoordinates(escr
452      if (!numSamplesEqual(&x, 1, getNumNodes()))      if (!numSamplesEqual(&x, 1, getNumNodes()))
453          throw RipleyException("setToX: Illegal number of samples in Data object");          throw RipleyException("setToX: Illegal number of samples in Data object");
454    
455        pair<double,double> xdx = getFirstCoordAndSpacing(0);
456        pair<double,double> ydy = getFirstCoordAndSpacing(1);
457      arg.requireWrite();      arg.requireWrite();
458  #pragma omp parallel for  #pragma omp parallel for
459      for (dim_t i1 = 0; i1 < m_N1; i1++) {      for (dim_t i1 = 0; i1 < m_N1; i1++) {
460          for (dim_t i0 = 0; i0 < m_N0; i0++) {          for (dim_t i0 = 0; i0 < m_N0; i0++) {
461              double* point = arg.getSampleDataRW(i0+m_N0*i1);              double* point = arg.getSampleDataRW(i0+m_N0*i1);
462              point[0] = (m_l0*(i0+m_offset0))/m_gNE0;              point[0] = xdx.first+i0*xdx.second;
463              point[1] = (m_l1*(i1+m_offset1))/m_gNE1;              point[1] = ydy.first+i1*ydy.second;
464          }          }
465      }      }
466  }  }
# Line 337  void Rectangle::assembleCoordinates(escr Line 468  void Rectangle::assembleCoordinates(escr
468  //private  //private
469  void Rectangle::populateSampleIds()  void Rectangle::populateSampleIds()
470  {  {
471      const index_t firstId = getNumNodes()*m_mpiInfo->rank;      // identifiers are ordered from left to right, bottom to top on each rank,
472      const index_t diff0 = m_N0*(m_N1-1)+1;      // except for the shared nodes which are owned by the rank below / to the
473      const index_t diff1 = m_N0*((m_NX-1)*m_N1+1);      // left of the current rank
474    
475        // build node distribution vector first.
476        // m_nodeDistribution[i] is the first node id on rank i, that is
477        // rank i owns m_nodeDistribution[i+1]-nodeDistribution[i] nodes
478        m_nodeDistribution.assign(m_mpiInfo->size+1, 0);
479        m_nodeDistribution[1]=getNumNodes();
480        for (dim_t k=1; k<m_mpiInfo->size-1; k++) {
481            const index_t x=k%m_NX;
482            const index_t y=k/m_NX;
483            index_t numNodes=getNumNodes();
484            if (x>0)
485                numNodes-=m_N1;
486            if (y>0)
487                numNodes-=m_N0;
488            if (x>0 && y>0)
489                numNodes++; // subtracted corner twice -> fix that
490            m_nodeDistribution[k+1]=m_nodeDistribution[k]+numNodes;
491        }
492        m_nodeDistribution[m_mpiInfo->size]=getNumDataPointsGlobal();
493    
494      m_nodeId.resize(getNumNodes());      m_nodeId.resize(getNumNodes());
495    
496        // the bottom row and left column are not owned by this rank so the
497        // identifiers need to be computed accordingly
498        const index_t left = (m_offset0==0 ? 0 : 1);
499        const index_t bottom = (m_offset1==0 ? 0 : 1);
500        if (left>0) {
501            const int neighbour=m_mpiInfo->rank-1;
502            const index_t leftN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
503    #pragma omp parallel for
504            for (dim_t i1=bottom; i1<m_N1; i1++) {
505                m_nodeId[i1*m_N0]=m_nodeDistribution[neighbour]
506                    + (i1-bottom+1)*leftN0-1;
507            }
508        }
509        if (bottom>0) {
510            const int neighbour=m_mpiInfo->rank-m_NX;
511            const index_t bottomN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
512            const index_t bottomN1=(neighbour/m_NX == 0 ? m_N1 : m_N1-1);
513  #pragma omp parallel for  #pragma omp parallel for
514            for (dim_t i0=left; i0<m_N0; i0++) {
515                m_nodeId[i0]=m_nodeDistribution[neighbour]
516                    + (bottomN1-1)*bottomN0 + i0 - left;
517            }
518        }
519        if (left>0 && bottom>0) {
520            const int neighbour=m_mpiInfo->rank-m_NX-1;
521            const index_t bottomN0=(neighbour%m_NX == 0 ? m_N0 : m_N0-1);
522            const index_t bottomN1=(neighbour/m_NX == 0 ? m_N1 : m_N1-1);
523            m_nodeId[0]=m_nodeDistribution[neighbour]+bottomN1*bottomN0-1;
524        }
525    
526        // the rest of the id's are contiguous
527        const index_t firstId=m_nodeDistribution[m_mpiInfo->rank];
528    #pragma omp parallel for
529        for (dim_t i1=bottom; i1<m_N1; i1++) {
530            for (dim_t i0=left; i0<m_N0; i0++) {
531                m_nodeId[i0+i1*m_N0] = firstId+i0-left+(i1-bottom)*(m_N0-left);
532            }
533        }
534    
535        for (dim_t k=0; k<m_nodeDistribution.size(); k++) {
536            cout << m_nodeDistribution[k] << endl;
537        }
538        cout<<endl;
539      for (dim_t k=0; k<getNumNodes(); k++) {      for (dim_t k=0; k<getNumNodes(); k++) {
540          index_t id = firstId+k;          cout << m_nodeId[k] << endl;
541          if (m_offset0 > 0 && k%m_N0==0)      }
542              id -= diff0;  
543          if (m_offset1 > 0 && k/m_N0==0)      // elements
544              id -= diff1;      m_elementId.resize(getNumElements());
545          m_nodeId[k]=id;  #pragma omp parallel for
546        for (dim_t k=0; k<getNumElements(); k++) {
547            m_elementId[k]=k;
548        }
549    
550        // face elements
551        m_faceId.resize(getNumFaceElements());
552    #pragma omp parallel for
553        for (dim_t k=0; k<getNumFaceElements(); k++) {
554            m_faceId[k]=k;
555      }      }
556  }  }
557    

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

  ViewVC Help
Powered by ViewVC 1.1.26